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Abstract 

The  purpose  of  this  research  effort  is  to  improve  and  characterize  range  estimation 
in  a  three-dimensional  FLASH  LAser  Detection  And  Ranging  (3D  FLASH  LADAR)  by 
investigating  spatial  dimension  blurring  effects.  The  myriad  of  emerging  applications  for 
3D  FLASH  LADAR  both  as  primary  and  supplemental  sensor  necessitate  superior  per¬ 
formance  including  accurate  range  estimates.  Along  with  range  information,  this  sensor 
also  provides  an  imaging  or  laser  vision  capability.  Consequently,  accurate  range  estimates 
would  also  greatly  aid  in  image  quality  of  a  target  or  remote  scene  under  interrogation. 

Unlike  previous  efforts,  this  research  accounts  for  pixel  coupling  by  defining  the 
range  image  mathematical  model  as  a  2D  convolution  between  the  system  spatial  impulse 
response  and  the  object  (target  or  remote  scene)  at  a  particular  range  slice.  Using  this  model, 
improved  range  estimation  is  possible  by  object  restoration  from  the  data  observations. 
Object  estimation  is  principally  performed  by  deriving  a  blind  deconvolution  Generalized 
Expectation  Maximization  (GEM)  algorithm  with  the  range  determined  from  the  estimated 
object  by  a  normalized  correlation  method.  Theoretical  derivations  and  simulation  results 
are  verified  with  experimental  data  of  a  bar  target  taken  from  a  3D  FLASH  LADAR  system 
in  a  laboratory  environment.  Simulation  examples  show  that  the  GEM  object  restoration 
improves  range  estimation  over  the  unprocessed  data  and  a  Wiener  filter  method  by  75% 
and  26%  respectively.  In  the  laboratory  experiment,  the  GEM  object  restoration  improves 
range  estimation  by  34%  and  18%  over  the  unprocessed  data  and  Wiener  filter  method 
respectively. 

This  research  also  derives  the  Cramer-Rao  bound  (CRB)  on  range  separation  estima¬ 
tion  of  two  point  sources  interrogated  by  a  3D  FLASH  LADAR  system.  Using  an  unbiased 
estimator,  range  separation  estimation  variance  was  attained  through  simulation  and  com¬ 
pared  favorably  to  the  range  separation  CRB  theory.  The  results  show  that  the  CRB  does 
indeed  provide  a  lower  bound  on  the  range  separation  estimation  variance  and  that  the  es- 
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timator  is  nearly  efficient.  With  respect  to  the  estimator,  traditional  pixel-based  estimators 
like  peak  detection  and  matched  filtering  are  biased  because  they  assume  there  is  only  one 
target  in  the  pixel.  Therefore,  an  unbiased  estimator  was  derived  accounting  for  the  possi¬ 
bility  of  two  targets  within  a  single  pixel. 

Additionally,  among  other  factors,  the  range  separation  CRB  is  a  function  of  two 
LADAR  design  parameters  (range  sampling  interval  and  transmitted  pulse-width),  which 
can  be  optimized  using  the  expected  range  resolution  between  two  point  sources.  Typi¬ 
cally,  a  shorter  transmitted  pulse-width  corresponds  to  better  range  resolution  (the  ability 
to  resolve  two  targets  in  range).  Given  a  particular  range  sampling  capability  determined 
by  the  receiver  electronics,  the  CRB  theory  shows  there  is  an  optimal  pulse-width  where 
a  shorter  pulse-width  would  increase  estimation  variance  due  to  the  under-sampling  of  the 
pulse  and  a  longer  pulse-width  would  degrade  the  resolving  capability.  Using  both  CRB 
theory  and  simulation  results,  an  investigation  is  accomplished  that  finds  the  optimal  pulse- 
width  for  several  range  sampling  scenarios.  For  example,  given  a  Gaussian  received  pulse 
model  sampled  every  1.876  ns,  both  range  separation  CRB  theory  and  simulation  predict 
an  optimal  pulse-width  standard  deviation  equal  to  0.88  ns.  As  the  speed  of  the  optical  re¬ 
ceiver  is  increased,  the  range  resolution  is  improved  with  a  corresponding  narrower  optimal 
pulse- width  attained  by  the  ability  to  sufficiently  sample  a  shorter  pulse-width.  Conversely, 
the  optimal  pulse- width  is  wider  with  slower  electronics  with  an  associated  negative  impact 
on  range  resolution. 
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Improving  Range  Estimation  of  a 


3-Dimensional  FLASH  LADAR 
via  Blind  Deconvolution 

I.  Introduction 

The  ability  to  accurately  view  a  remote  scene  has  long  been  a  human  military  en¬ 
deavor.  From  primeval  warriors  using  mountains  or  trees  to  see  troop  formations 
with  their  naked  eye  to  early  seafarers  using  primitive  optics  to  assess  ship  capabilities  or 
harbor  defenses  to  today’s  combatants  operating  advanced  optics  (manned  and  unmanned 
platforms)  and  imaging  satellites  to  observe  troop  or  missile  movements,  the  advantage  to 
the  military  that  can  accurately  assess  the  remote  battlefield  has  never  been  questioned. 

With  modem  technology  development,  remote  sensing  has  advanced  and,  in  one  par¬ 
ticular  sensor  area,  has  bonded  with  another  indispensable  military  capability:  RAdio  De¬ 
tection  And  Ranging  (RADAR).  Since  World  War  II,  RADAR  capability  has  been  a  critical 
technology  with  respect  to  offensive  and  defensive  capabilities  and  missile  defense.  Ad¬ 
vances  in  RADAR  have  steadily  progressed  since  the  early  effective  use  of  RADAR  by 
Great  Britain  against  Germany  in  the  Battle  of  Britain.  However,  RADAR  is  fundamentally 
limited  in  some  ways  by  its  operating  wavelength  in  the  electromagnetic  spectrum.  One 
of  the  latest  advancements  is  in  the  field  of  RADAR  is  adapting  the  use  of  lasers  to  the 
ranging  issue  and  developing  a  LAser  Detection  And  Ranging  (LADAR)  system.  LADAR 
allows  for  the  benefit  of  a  smaller  operating  wavelength  (e.g.,  resolution  and  material  inter¬ 
action)  and  the  directionality  of  laser  transmissions  (which  are  ideal  for  urban  environment 
interrogation)  while  still  retaining  the  imaging  and  ranging  benefit  of  a  traditional  RADAR. 
Just  as  there  is  no  one  branch  of  the  military  that  can  operate  independently  in  the  mod¬ 
ern  battlefield,  a  LADAR  is  not  meant  to  be  a  panacea  and  there  are  applications  where 
RADAR  is  still  preferred.  Though,  as  the  technology  continues  to  mature,  LADAR  will  be 
an  invaluable  contributor  in  imaging  and  ranging  sensor  suite  available  to  the  warfighter. 


1 


Motivation  for  this  research  effort  is  accomplished  in  this  chapter  by  introducing 
LADAR  and  explaining  the  importance  of  range  estimation.  Specifics  of  the  research  con¬ 
tributions  are  described  with  corresponding  benefits.  Finally,  the  document  organization  is 
given. 

1.1  Motivation 

The  driving  force  behind  this  research  endeavor  is  to  adopt  a  realistic  physical  model 
for  the  return  signals  of  a  three-dimensional  FLASH  LAser  Detection  And  Ranging  (3D 
FLASH  LADAR)  and  then  develop  methods  of  improving  the  most  vital  unknown  param¬ 
eter  from  that  model:  range  to  target.  More  precise  range  measurements  aid  intelligence 
gathering,  target  recognition,  mapping,  imaging,  object  classification,  navigation,  and  pre¬ 
cision  strike  capabilities.  The  trend  towards  computer  vision  systems  with  active  illumina¬ 
tion  necessitates  the  use  of  3D  FLASH  LADARs  capable  of  rapid  range  data  acquisition 
with  a  wide  enough  field  of  view  (FOV)  to  allow  the  operator  access  to  an  appropriate  bat¬ 
tlefield  representation.  By  acquiring  the  remote  scene  in  this  manner,  however,  the  sensor 
will  be  negatively  affected  by  the  spatial  blurring  inherent  in  the  image  formation  process. 

The  importance  of  being  able  to  correctly  range  to  the  remote  environment  is  charac¬ 
terized  by  General  T.  Michael  Moseley  in  a  2007  CSAF  white  paper: 

The  Air  Force  is  often  first  to  the  fight  and  last  to  leave.  We  give  unique  options 
to  all  Joint  Force  Commanders.  The  Air  Force  must  safeguard  our  ability  to: 
see  anything  on  the  face  of  the  earth;  range  it;  observe  or  hold  it  at  risk;  supply, 
rescue,  support  or  destroy  it;  assess  the  effects;  and  exercise  global  command 
and  control  of  all  these  activities.  Rising  to  the  21st  Century  challenge  is  not  a 
choice.  It  is  our  responsibility  to  bequeath  a  dominant  Air  Force  to  Americas 
joint  team  that  will  follow  us  in  service  to  the  Nation  [57]. 

Unlike  3D  scanning  LADARs  that  build  a  3D  scene  by  rastering  multiple  laser  scans 
with  a  dwell  required  at  each  point,  a  3D  FLASH  LADAR  system  produces  a  set  of  se¬ 
quential  two-dimensional  (2D)  images  due  to  a  fast  range  gate  (i.e.  shutter)  resulting  in  a 
three-dimensional  data  cube  (spatial  and  range)  of  the  remote  scene.  In  reality,  the  sensor 
captures  a  fourth  dimension  which  is  the  photo-electron  count  for  each  volume  element 
(voxel).  Each  2D  range  slice  image  contains  the  detected  photo-electrons  at  each  pixel  for 
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a  particular  range.  The  photo-electron  counts  are  directly  proportional  to  the  return  signal 
intensities  incident  upon  the  detector.  Unique  to  the  FLASH  LADAR  sensors,  each  pixel 
in  the  array  detects  its  own  attenuated  and  time-delayed  version  of  the  transmitted  signal. 
Investigating  pixel  data,  the  blurring  effects  are  evident  in  the  pixel’s  received  waveform. 

3D  FLASH  LADAR  range  estimation  errors  of  a  remote  scene  can  occur  due  to  sev¬ 
eral  system  factors  including  the  optical  spatial  impulse  response  (diffraction  limited,  atmo¬ 
spheric  turbulence),  detector  blurring,  photon  noise,  and  readout  noise.  These  factors  either 
cause  the  scene’s  intensity  to  spread,  or  blur,  across  pixels  or  add  unwanted  and  disruptive 
noise  effects.  The  intensity  spreading  and  noise  corrupts  the  correct  pixel  intensities  at  each 
range  gate  by  mixing  intensities  with  neighboring  pixels  thereby  providing  false  intensity 
values  and  therefore  incorrect  photon  counts  to  the  range  estimator.  Without  blur  and  noise 
compensation,  the  range  estimates  would  then  be  inaccurate  to  a  degree  depending  on  the 
blur  and  noise  severity. 

3D  FLASH  LADAR’s  popularity  is  increasing  due  to  its  small  size,  rapid  image  ac¬ 
quisition,  and  range  resolving  capabilities.  There  are  several  examples  that  highlight  appli¬ 
cations  in  practical  situations:  As  part  of  its  return  to  flight  efforts  following  the  Columbia 
disaster,  NASA  uses  a  3D  imaging  LADAR  to  inspect  the  integrity  of  the  Space  Shuttle’s 
Thermal  Protection  System  prior  to  reentry  [45].  Sandia  National  Laboratory  developed  a 
counter-sniper  3D  coherent  detection  LADAR  sensor  designed  to  trace  the  source  of  the 
bullet  by  optical  signature  and  bullet  trajectory  analysis  [74]. 

Augmenting  the  FLASH  technology  to  the  LADAR’s  active  sensing  capability,  the 
possibilities  of  future  technology  include  remote  video  feeds  from  airborne  or  spacebome 
platforms  to  command  and  control  centers,  precise  autonomous  navigation  in  GPS-denied 
regions,  autonomous  precision  strike  with  guided  cruise  missiles  or  intelligent  gravity  mu¬ 
nitions,  and  battlefield  awareness  in  day/night  conditions  for  airborne  or  ground  forces  in 
dynamic  environments. 
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1.2  3D  FLASH  LADAR  Research  Contributions 


1.2.1  Improving  Range  Estimation  by  Spatial  Processing  ( Chapter  V)  .  Previ¬ 
ous  work  in  3D  FLASH  LADAR  has  only  modeled  an  ideal  return  per  pixel  and  not  the 
real  world  effects  of  spatial  blurring  [9],  [38].  This  research  will  enhance  the  model  by 
adding  the  spatial  impulse  response  thereby  considering  all  the  pixel’s  signals  in  the  range 
estimation  algorithm  for  a  particular  pixel.  The  benefit  of  this  research  is  for  future  imple¬ 
mentation  in  an  operational  environment.  Previously,  a  3D  representation  of  a  remote  scene 
was  built  by  single-pixel  LADAR  scanners.  Consequently,  the  scanning  3D  LADARs  have 
limited  spatial  extents  on  each  collect  and  do  not  see  the  effects  of  spatial  blurring.  As 
laser  vision  hardware  improves,  the  development  will  trend  towards  FLASH  systems  that 
capture  scene  data  very  rapidly  over  a  large  pixel  array.  Given  proper  spatial  sampling,  this 
method  of  data  capture  would  see  the  effects  of  spatial  blurring.  The  spatial  blurring  would 
contribute  negatively  to  current  methods  of  range  estimation  because  each  pixel’s  return 
waveform  would  interact  with  those  of  its  neighboring  pixels.  New  estimation  solutions 
must  be  developed  that  account  for  these  blurring  effects  and  essentially  “deblur”  the  data 
to  increase  range  estimation  performance.  This  research  builds  this  enhanced  model  and 
improves  range  estimation  by  spatially  processing  the  data  using  a  well-known  spatial  filter 
(deconvolution)  and  a  novel  object  recovery  algorithm  (blind  deconvolution)  [56]. 

1.2.2  Unbiased  Two  Point  Target  Temporal  and  Spatial  Estimator  (Chapter  IV). 
This  contribution  supports  the  CRB  work  from  the  previous  section.  Given  the  two  target 
model,  conventional  pixel-based  estimators  like  peak  detection  and  matched  filtering  are 
biased  because  they  assume  there  is  only  target  in  the  pixel.  Therefore,  an  unbiased  es¬ 
timator  was  developed  accounting  for  the  possibility  of  two  targets  within  a  single  pixel. 
Based  on  a  least  sum  squares  approach,  the  ability  to  sufficiently  estimate  the  ranges  and 
amplitudes  of  two  point  targets  is  developed  and  verified  to  be  unbiased. 

1.2.3  Lower  Bound  on  Range  Separation  Estimation  Variance  ( Chapter  VI)  .  An 
important  metric  of  a  physical  model  with  several  unknowns  is  to  understand  the  optimal  es- 
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timator  variance  achievable  regardless  of  the  specific  estimation  scheme.  The  Cramer-Rao 
Bound  (CRB)  provides  the  lower  bound  on  estimator  variance  given  an  unbiased  estimator. 
Previous  CRB  work  in  3D  FLASH  LADAR  adopted  a  physical  model  that  did  not  account 
for  the  spatial  blurring  between  pixels  [9],  [38].  The  benefit  of  including  these  spatial 
effects  in  the  CRB  development  is  that  the  estimation  and  CRB  results  would  now  be  nega¬ 
tively  affected  by  the  signals  from  adjacent  pixels  to  a  degree  depending  on  the  pixel  range 
differences.  A  two  point  target  scene  model  is  adopted  to  show  the  CRB  on  range  separa¬ 
tion  estimation.  The  effects  of  changing  the  separation  are  shown  to  drastically  affect  the 
ability  to  estimate  that  separation. 

1.2.4  Optimal  Pulse-Width  based  on  Range  Resolution  (Chapter  VI)  .  Utilizing 
the  CRB  and  unbiased  two  point  target  range  separation  estimator,  a  method  is  developed 
where  an  optimal  pulse-width  is  determined  based  on  the  expected  range  resolution  using 
the  two  point  target  model.  Typically,  a  shorter  transmitted  pulse-width  corresponds  to  bet¬ 
ter  range  resolution  (the  ability  to  estimate  two  distinct  targets  in  range).  Given  a  particular 
range  sampling  capability  determined  by  the  receiver  electronics,  the  CRB  and  simulation 
shows  there  is  an  optimal  pulse-width  where  a  shorter  pulse-width  would  increase  esti¬ 
mation  variance  due  to  the  under-sampling  of  the  pulse  and  a  longer  pulse-width  would 
degrade  the  resolving  capability.  Using  two  distinct  and  separate  techniques  of  CRB  and 
simulation,  an  investigation  is  accomplished  that  finds  the  optimal  pulse-width  for  several 
range  sampling  scenarios.  Benefits  of  this  analysis  include  the  ability  to  aid  in  LADAR 
system  design  using  independent  statistical  methods  (CRB). 

1.3  Organization 

The  dissertation  is  organized  as  follows:  Chapter  II  provides  background  theory,  data 
model,  and  a  discussion  of  previous  LADAR  research.  Chapter  III  details  the  3D  FLASH 
LADAR  hardware  used  in  experimental  collects  as  well  as  the  procedures  used  to  condition 
the  data  for  appropriate  use  for  the  selected  mathematical  model.  Chapter  IV  contains  the 
pertinent  pixel-based  range  estimation  algorithms.  Chapter  V  shows  that  object  recovery 
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does  improve  range  estimation.  Chapter  VI  derives  the  CRB  for  range  separation  estima¬ 
tion  and  predicts  an  optimal  pulse-width  that  provides  the  best  range  resolution.  Finally, 
Chapter  VII  summarizes  the  research  contributions  and  outlines  future  research  opportuni¬ 
ties. 
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II.  Background 


This  chapter  serves  as  a  review  of  background  theory  and  previous  research  related 
to  three-dimensional  FLASH  LAser  Detection  And  Ranging  (3D  FLASH  LADAR). 
The  focus  on  the  theory  and  literature  review  will  be  related  to  the  major  topic  areas:  range 
estimation,  spatial  processing,  performance  bounding,  and  optimal  parameter  selection. 

This  chapter  is  organized  as  follows:  Section  2.1  discusses  imaging  and  coherence 
theory  and  how  it  applies  to  3D  FLASH  LADAR.  Sections  2.2,  2.3,  and  2.4  discuss  de- 
convolution,  maximum  likelihood  parameter  estimation,  and  the  Generalized  Expectation 
Maximization  (GEM)  algorithm  respectively.  Section  2.5  describes  the  data  model  that  will 
be  used  in  Chapters  III,  IV,  and  V.  Finally,  Section  2.6  reviews  previous  research  related 
to  LADAR  data  processing,  blind  deconvolution,  bounding  performance,  and  parameter 
optimization. 

2.1  LADAR  Imaging  Theory 

The  goal  of  this  section  is  to  describe  the  3D  FLASH  LADAR  imaging  operation 
as  a  linear  and  spatially-invariant  system.  Linear  systems  theory  has  many  benefits  with 
the  chief  benefit  of  being  able  to  describe  the  observed  data  (image)  as  a  convolution  of 
the  object’s  intensity  with  a  spatial  impulse  response.  This  convolution  is  an  integral  part 
of  the  mathematical  model  used  in  this  research  describing  the  detected  photons  in  the  3D 
FLASH  LADAR.  The  spatial  impulse  response  completely  describes  the  optical  system  to 
include  any  random  atmospheric  disturbance.  The  argument  that  optical  imaging  can  be 
cast  in  the  linear  system  framework  has  been  established  in  the  literature  [24],  [25].  Similar 
arguments  are  made  here  to  verify  that  this  framework  is  applicable  to  this  research.  A 
foundational  understanding  of  why  this  object-image  relationship  holds  is  key  because  it 
allows  the  use  of  object  reconstruction  algorithms  from  the  simple  inverse  filter  to  the  more 
complicated  blind  deconvolution  methods.  First,  a  method  is  needed  to  accurately  describe 
the  illuminating  light’s  movement  and  interaction  with  its  environment  and  how  light  prop¬ 
agates.  Following,  the  linear  system  framework  can  be  constructed  with  an  example  of  a 
spatial  impulse  response  for  a  simple  imaging  system. 
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2.1.1  Description  of  Light.  From  [30],  there  are  three  mathematical  descriptions 
in  which  the  light  used  in  optical  imaging  systems  can  be  described.  They  are  geometrical 
optics,  quantum  optics,  and  physical  (wave)  optics. 

The  simplest  and  least  accurate  mathematical  model  is  called  geometrical  optics  (GO) 
and  is  a  good  approximation  when  wavelengths  are  small  compared  to  the  dimensions  of 
the  optics.  GO  analysis  operates  on  the  principle  of  light  described  as  rays  and  is  a  valid 
technique  to  determine  basic  properties  of  an  imaging  system  like  object  distance,  magnifi¬ 
cation,  and  pixel  area  at  the  target.  For  instance,  the  location  of  an  object’s  image  through 
an  imaging  system  is  based  on  the  ray  tracing  method.  While  being  able  to  show  optical 
abberations,  GO  does  not  handle  diffraction  or  interference  effects  and  predicts  the  loca¬ 
tion  of  an  image  to  be  a  point  (without  aberrations).  The  blurring  of  the  imaging  system  is 
not  accounted  for  in  GO  which  makes  it  a  poor  choice  to  describe  light  propagation  in  this 
research. 

From  [30],  the  most  accurate  and  complex  mathematical  description  of  light  is  called 
quantum  optics  (QO)  and  is  valid  in  all  optical  scenarios  (wavelength,  irradiance  levels, 
and  optic  dimensions).  In  QO,  light  is  considered  an  electromagnetic  wave  with  its  en¬ 
ergy  quantized  into  massless  particles  called  photons  rather  than  a  continuous  wave.  While 
being  the  most  physically  accurate,  computations  tend  to  be  slow  and  cumbersome.  Con¬ 
cerning  imaging  applications,  the  extra  time  and  resources  required  for  QO  is  not  beneficial 
when  trying  to  understand  and  mitigate  the  macroscopic  light  blurring  effects  of  an  optical 
system. 

Also  from  [30],  the  remaining  mathematical  method  is  called  physical  optics  (PO), 
or  wave  optics,  in  which  the  light  is  considered  to  propagate  as  a  transverse  electromag¬ 
netic  wave.  In  general,  PO  can  be  used  to  describe  diffraction  and  interference  effects  by 
Maxwell’s  equations  using  a  scalar  theory  approximation.  Rather  than  a  vector-based  the¬ 
ory,  scalar  theory  is  valid  for  describing  light  behavior  when  the  wavelength  of  the  light 
is  much  less  than  the  dimensions  of  the  diffracting  objects  and  when  traveling  through  a 
uniformly  dielectric  medium.  Whereas  GO  describes  an  image  to  be  a  point,  PO  utilizes 
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Object  Plane  Focal  Ptane  Image  Plane 

Figure  2.1:  Simplified  depiction  of  an  imaging  system.  In  the  3D  FLASH  LADAR,  the 
object  is  the  target  under  illumination,  the  lens  are  the  front  end  optics,  and 
the  image  would  be  the  2D  range  slice  image  of  the  intensity  return  from  the 
object  at  a  particular  return  time. 

diffraction  effects  through  the  optical  system  to  depict  the  image  spread  about  the  point  that 
GO  predicted.  While  there  may  be  situations  where  the  physical  optics  assumptions  and 
approximations  fail,  the  more  accurate  quantum  optics  approach  tends  to  be  impractical 
due  to  increased  complexity  and  processing  times  for  operational  use.  In  most  practical 
situations,  PO  is  sufficient  to  describe  the  light’s  movement  and  interactions  with  structures 
given  high  enough  irradiance  levels. 

A  common  practice  in  imaging  systems  is  to  treat  light  as  an  electromagnetic  wave 
using  PO  until  the  light  hits  the  detector  in  which  the  light  is  then  considered  a  particle 
or  photon.  This  assumption  allows  for  the  benefit  of  the  dual  nature  of  light  and  will  be 
adopted  for  this  research.  Furthermore,  PO  is  sufficiently  accurate  to  describe  an  optical 
imaging  system  as  a  linear  system. 

2.1.2  Optical  Field  Propagation.  Based  on  [25]  and  referring  to  Figure  2.1,  the 
purpose  of  this  section  is  to  be  able  to  describe  how  to  mathematically  propagate  an  optical 
field  from  one  plane  to  another  with  varying  levels  of  accuracy.  In  order  to  mathematically 
propagate  an  optical  field,  a  diffraction  formula  must  be  used  due  to  the  many  point  sources 
in  the  observation  plane. 
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Through  Maxwell’s  equations,  the  Huygens-Fresnel  principle,  and  Kirchhoff’s  the¬ 
ory,  a  closed-form  mathematical  solution  for  the  optical  field  at  a  remote  point  can  be 
attained  called  the  Rayleigh-Sommerfeld  diffraction  formula.  This  diffraction  formula 
is  a  general  result  from  scalar  diffraction  theory  with  the  only  assumptions  lying  within 
scalar  diffraction  theory.  With  monochromatic  and  narrowband  assumptions,  the  Rayleigh- 
Sommerfeld  diffraction  formula  is  given  by  the  following  equation  for  the  complex  phasor 
Ua  of  the  scalar  optical  field  at  a  distance  z  away  from  the  source  field  U  [25]: 


Ua(u,v)  = 


TT/J.  s  j2nR(j,TI,U,v) 

U(t,v)e  A 

3  *R(£,V,u,v) 


d£drj 


— oo  — oo 


(2.1) 


where  (u,  v )  are  observation  plane  coordinates,  (£,  rf)  are  source  plane  coordinates,  A  is  the 
mean  wavelength,  and  i?(£,  77,  u,  v )  =  \J z2  +  (£  —  u )2  +  (77  —  v )2  is  the  distance  between 
every  point  in  the  source  plane  to  every  point  in  the  observation  plane.  The  complex  phasor 
is  related  to  the  scalar  optical  field  by 


ua  (u,  v,  t )  =  Re  {Ua  (u,  v )  exp  (-j2-Kist)}  (2.2) 

where  Re  {}  means  the  “real  part”,  j  =  \f—l,  v  is  the  frequency  of  the  light,  and  t  signifies 
time.  The  optical  field  theory  focuses  on  the  complex  phasor  U  development  since  the  time 
dependence  is  already  known  [25]. 

A  useful  simplification  called  the  Frensel  approximation  (near-field  or  paraxial  ap¬ 
proximation)  can  be  employed  to  reduce  the  complexity  of  the  range  term,  although  the 
two  instances  of  the  range  term  need  to  be  handled  differently.  Small  errors  in  range  term 
in  the  denominator  are  usually  not  critical  due  to  the  range  to  target  (z)  being  much,  much 
bigger  than  the  spatial  extents  in  the  observation  and  image  plane.  Conversely,  small  errors 
in  the  range  term  residing  in  the  exponential  can  be  significant  given  that  it  is  divided  by  the 
light’s  wavelength  which  is  on  the  order  of  hundreds  of  nanometers  in  the  light  or  infrared 
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spectrum.  Using  the  binomial  approximation  for  R  in  the  exponential  results  in 


R  =  \J ' z2  +  (£  -  uf  +  (rj  -  v f 

=  zyj 1  +  ((^  -  u)  /zf  +  ((r/  -  v)  /zf 
=  z{  1  +  0.5  ((£  -  «)  /zf  +  0.5  ((77  -  v)  /zf) 

~  „  ,  (£-«)2  ,  (v~vf 


(2.3) 


and  approximating  R  pz  z  in  the  denominator,  results  in  the  Fresnel  diffraction  formula 
given  by  [25] 

Ua(M,n)  =  ^-y  J  U(Z,rj)ej^-u)2+{r]-v)2}dZdr}  (2.4) 

—  OO  —OO 


which  describes  a  convolution  operation  for  the  free  space  propagation  of  an  optical  field 
from  one  plane  to  another.  A  spatial  impulse  response  (spatial  point  spread  function)  for 
free  space  propagation  is  then  defined  by 


j2nz 


(2.5) 


It  is  interesting  to  note  that  even  free-space  propagation  can  be  cast  in  the  linear  systems 
framework.  A  later  section  will  make  the  argument  that  an  imaging  system  can  be  repre¬ 
sented  as  a  linear  system  as  well. 

An  alternate  way  to  view  the  free-space  Fresnel  diffraction  integral  is  by  factoring 
out  the  variables  that  don’t  depend  on  the  variables  of  integration  and  results  in 


Ua(u,v )  = 


']2ttz 
£  A  £ 


jXz 


_T/  .  N  J7r(g2+T72)  -j2n(u£+vri) 

t/(£,  r])e  Xz  e  Xz  d£dr). 


—00  —00 


(2.6) 


which  is  a  scaled  Fourier  Transform  of  the  aperture  field  and  the  quadratic  exponential.  The 
Fresnel  diffraction  formula  still  accounts  for  the  curvature  of  the  wavefront,  but  assumes 
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a  parabolic  rather  than  spherical  wavefront  shape.  While  Equation  (2.6)  specifies  the  op¬ 
tical  field  (volts/meter)  at  a  distance,  the  intensity  at  that  point  is  the  quantity  of  interest 
in  imaging.  Considering  the  wave  is  monochromatic,  the  intensity  (watts/meter2)  can  be 
determined  by  taking  the  magnitude  squared  of  the  complex  phasor  of  the  optical  field  or 
Ia  (u,  v )  =  \Ua  {u,  v )  | 2 .  When  the  wave  is  not  monochromatic,  the  intensity  becomes  the 
time-average  (•)  of  the  amplitude  squared  of  the  scalar  optical  field 

Ia(u,v)  =  (\ua(u,v,t)\2)  (2.7) 

where  ua  (?./,  v.  t )  was  defined  in  Equation  (2.2).  All  future  references  to  an  “optical  field” 
refer  to  the  complex  phasor  U  unless  explicitly  stated  otherwise. 

Equation  (2.4)  or  (2.6)  can  now  be  used  to  describe  the  imaging  operation  where 
optics  are  placed  between  the  object  and  image.  The  next  section  summarizes  the  resulting 
impulse  response  of  a  general  imaging  scenario. 

2.1.3  Impulse  Response  of  an  Imaging  System  with  a  Thin  Lens.  The  purpose  of 
this  subsection  is  to  illustrate  an  example  of  the  impulse  response  from  a  simple  imaging 
architecture.  The  imaging  system  converts  the  diverging  spherical  waves  emanating  from 
an  object  to  converging  spherical  waves  culminating  at  the  image.  The  lens  is  assumed 
to  be  a  thin  lens  meaning  the  light  enters  and  leaves  the  lens  at  the  same  coordinates.  Of 
course,  there  is  a  diameter  to  all  lenses  as  well  as  irregularities  that  make  this  assumption 
invalid.  However,  it  will  suffice  for  the  purposes  of  a  theoretical  understanding  of  the  lens’ 
effect  on  incident  light. 

In  general,  the  purpose  of  the  imaging  system  is  to  reproduce  an  object  in  a  better 
manner  than  possible  without  the  system.  With  no  aberrations,  the  geometrical  optics  anal¬ 
ysis  predicts  a  “perfect”  image  aside  from  a  scaling  term,  although  this  image  is  only  valid 
as  the  wavelength  goes  to  infinity  (A  — >•  0).  Wave  optics  predicts  a  more  physically  accurate 
image  that  is  dominated  by  the  effects  of  diffraction.  As  stated  previously,  a  significant  con¬ 
cept  in  this  research  is  that  the  3D  LADAR  is  operating  in  a  linear  system.  This  assumption 
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allows  for  the  LADAR  to  be  entirely  represented  by  a  spatial  impulse  response.  The  images 
are  then  produced  by  the  convolution  of  the  object  and  spatial  impulse  response.  The  key  is 
to  be  able  to  describe  an  optical  imaging  system  by  a  spatial  impulse  response.  By  placing 
a  point  source  in  front  of  a  lens,  the  impulse  response  of  the  lens  can  be  attained.  This 
lens  impulse  response  is  valid  for  compound  or  more  complex  optics  since  all  the  imaging 
system  optics  convert  a  diverging  spherical  wave  into  a  converging  spherical  wave. 

Under  the  general  assumption  of  the  linearity  of  wave  propagation,  the  relationship 
between  a  field  at  the  image  and  object  plane  can  be  given  by  a  superposition  integral  [25]: 

OO  OO 

Ui  (u, v)  —  J  I  h  (u,  v;  £,  77)  U0  (£,  rj)  d^drj  (2.8) 

— 00  — OO 

where  (Jt  and  UQ  are  the  image  and  object  plane  optical  field  respectively  and  h  is  the 
impulse  response  and  is  an  optical  field  at  (u,  v )  produced  by  an  amplitude  point  source 
at  (£,77).  The  spatial  impulse  response  can  describe  optical  systems  from  simple  free- 
space  to  the  most  complicated  optics.  If  the  system  is  considered  space  invariant  (i.e.  an 
isoplanatic  imaging  situation  exists)  then  h  is  h  (u  —  £,  v  —  rf)  where  Equation  (2.8)  is  now 
a  convolution  integral.  From  [25],  however,  the  Fresnel  diffraction  integral  (Equation  (2.4)) 
is  used  along  with  the  phase  transformation  of  a  lens  to  derive  the  general  form  for  a  spatial 
impulse  response  of  a  single  thin  lens  to  be 


h  (u,  v\ £,  rj)  = 


1  i^^w)e^-{ew) 


e  z2 
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(2.9) 


the  wavelength,  z\  as  the  distance  from  the  object  to  focal  plane,  z2  as  the  distance  from  the 
focal  to  image  plane,  /  being  the  focal  length,  P  as  the  pupil  function,  and  (x,  y )  as  the  focal 
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plane  coordinates.  Using  assumptions  about  the  quadratic  phase  terms  and  normalizing  the 
coordinates  to  eliminate  effects  of  inversion  and  magnification,  the  general  form  reduces  to 
a  normalized  point  spread  function 


h(u,v)  =  - —$[P(x,y)\ruJij=v 


A  z- 


y—\zi 


(2.10) 


with  A  being  the  optical  field  amplitude,  S'  as  the  Fourier  transform  operator,  and  fy) 
are  the  focal  plane  spatial  frequency  coordinates.  It  is  only  under  specific  conditions  that 
Equation  (2.10)  results  from  the  more  general  impulse  response.  First,  the  lens  law  must  be 
satisfied: 


1  _  1  1 

7  A  z2 


(2.11) 


which  is  a  mandatory  condition  for  imaging  to  occur  and  thus 


(2.12) 


in  Equation  (2.9)  reduces  to  unity.  Second,  since  the  goal  of  imaging  is  to  obtain  the  in¬ 
tensity  of  the  image,  any  multiplicative  phase  terms  with  dependence  only  on  image  plane 
coordinates  can  be  discarded.  In  other  words,  the  term  exp  (u2  +  v2)  j  can  be  ignored. 
Finally,  the  quadratic  phase  term  dependent  on  object  plane  coordinates,  exp  (£2  +  r/2)  j , 
is  ignored  by  noting  that  the  object  is  a  point  source  and  the  span  of  the  object  coordinates 
are  very  small.  Therefore,  it  would  contribute  a  trivial  amount  to  the  intensity  on  the  focal 
plane.  With  these  three  conditions  satisfied,  the  impulse  response  for  a  thin  lens  takes  the 
form  of  Equation  (2. 10). 

This  result  is  an  example  of  an  ideal  impulse  response  for  an  optical  imaging  system. 

It  is  ideal  in  the  sense  that  there  are  no  aberrations  or  atmospheric  turbulence.  Using  the 
principle  planes  concept  from  geometrical  optics,  most  optics  in  an  imaging  systems  can 
be  considered  a  “thin  lens”  with  light  entering  the  system  with  one  orientation  and  exiting 
at  another  orientation  without  regard  to  the  inner  optical  structures.  Thus,  the  thin  lens 
impulse  response  is  a  good  approximation  or  starting  point  for  reconstruction  algorithms. 
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2.1.4  Optical  Imaging  as  a  Linear  and  Nonlinear  System.  In  order  to  validate 
the  mathematical  model  adopted  in  this  research,  the  relationship  between  the  object  in¬ 
tensity  and  image  intensity  needs  to  be  a  linear  relationship.  Depending  on  the  coherence 
properties  of  the  illuminating  light,  this  linear  relationship  may  or  may  not  exist.  This  sub¬ 
section  gives  examples  of  both.  The  next  subsection  concentrates  on  coherence  theory,  how 
it  affects  the  object-image  linear  relationship,  and  why  this  research  can  assume  a  linear 
relationship  between  the  object  and  image  intensity  does  exist  in  a  3D  FLASH  LADAR. 

Presented  again  for  convenience,  the  relationship  between  a  field  at  the  image  and 
object  plane  can  be  given  by  a  superposition  integral  due  to  the  linearity  of  wave  propaga¬ 
tion 

CO  CO 

Ui(u,v)=J  I  h  (u,  v;  £,  77)  U0  (£,  rj)  d^dr]  (2.13) 

—  CO  —CO 

where  Ui  and  U r>  are  the  image  and  object  plane  optical  fields  respectively  and  h  is  the 
impulse  response  and  is  an  optical  field  at  (u,  v)  produced  by  an  amplitude  point  source 
at  (£,77).  Again,  if  the  system  is  considered  space  invariant  (i.e.  an  isoplanatic  imaging 
situation  exists)  then  h  is  h  (u  —  £,v  —  rj)  where  Equation  (2.13)  is  now  a  convolution 
integral. 

In  a  simplified  imaging  situation,  the  imaging  system  consists  of  an  object,  a  lens, 
and  an  image.  The  ideal  image  predicted  by  geometrical  optics  is 


U9(u,v) 


U  V  \ 

M’  M/ 


(2.14) 


where  M  is  the  magnification  and  Ua  is  the  object.  This  ideal  image  is  the  result  of  the 
superposition  integral  as  A  — >  0.  Using  this  result  as  the  object  plane  amplitude  in  Equa¬ 
tion  (2.13),  the  field  at  image  plane  is  a  convolution  of  the  impulse  response  and  image 
predicted  by  geometrical  optics  [25]: 


Ui(u,  v )  =  h(u,  v )  <g)  Ug(u,  v ) 


(2.15) 
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This  result  highlights  the  spreading  effect  that  diffraction  imposes  on  the  ideal  image. 

Unless  further  propagation  is  necessary  where  the  optical  field  is  required,  the  optical 
intensity  at  the  detector  is  the  quantity  of  interest.  Results  for  image  intensity  will  be  stated 
here  and  justified  in  the  next  section  using  coherence  theory.  The  image  intensity  is  the 
time-averaged,  magnitude  squared  of  the  field  and  is  defined  by  an  intensity  convolution 
for  incoherent  illumination 


Ii(u,v)  =  \h(u,v)\2  ®  \Ug(u,v)\2  .  (2.16) 

This  result  for  image  intensity  is  the  important  result  of  this  section.  It  serves  as  the  ba¬ 
sis  for  the  mathematical  model  and  allows  for  advanced  techniques  for  object  restoration. 
If  coherent  illumination  is  encountered,  the  linear  relationship  for  intensity  vanishes  and 
image  intensity  is  defined  by  an  amplitude  convolution: 

Ii(u,v )  =  | h(u,v)  (8)  Ug(u,v) |2  .  (2.17) 

which  results  in  a  non-linear  relationship  between  the  object  and  image.  Clearly,  it  is  seen 
that  incoherent  illumination  is  linear  in  intensity  and  coherent  illumination  is  linear  in  am¬ 
plitude.  The  spatial  impulse  response  for  incoherent  illumination  is  the  amplitude  squared 
of  the  coherent  illumination  spatial  impulse  response. 

It  must  be  shown  or  proven  that  the  3D  FLASH  LADAR  produces  or  approaches 
incoherent  object  illumination  in  order  to  develop  algorithms  for  the  recovery  of  the  original 
object,  Ug,  using  deconvolution  algorithms.  Otherwise,  the  mathematical  model  would 
change  from  object  intensity  (i.e.  photon  counts)  to  object  field  recovery  in  order  to  benefit 
from  linear  systems  theory.  Since  the  observed  data  is  based  on  the  image  intensity,  backing 
out  the  object  field  from  coherent  illumination  would  require  other  methods  rather  than 
deconvolution. 

2.1.5  Coherence  Theory  and  Laser  Light  Statistics.  Using  [24]  and  [25],  this 
section  serves  as  background  on  coherence  theory  and  how  to  use  this  theory  to  express  the 
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image  intensity  as  in  Equation  (2.16)  or  (2.17).  Coherence  theory  also  dictates  the  statistics 
that  govern  the  laser  light  incident  on  the  detector  surface. 

The  image  intensity  related  to  different  types  of  coherence  is  governed  by  the  laser 
light’s  spatial  coherence  between  two  points  called  mutual  intensity.  In  order  to  understand 
how  coherence  affects  imaging,  the  monochromatic  light  assumption  has  to  be  relaxed  and 
the  light  model  changed  to  polychromatic.  This  yields  a  generic  optical  scalar  field  defined 
by 

u  (u,  v,  t )  =  {U  (u,  v ,  t)  exp  (—j2nut)}  .  (2.18) 


where  the  complex  phasor  is  changed  to  be  time-varying.  The  image  plane  complex  phasor 
Ui  (u,  v ,  t)  results  from  a  convolution  between  the  impulse  response  and  the  object  plane 
complex  phasor  Ug(£,r),t).  Neglecting  the  different  time  delays  from  different  coordinates, 
the  subsequent  image  plane  intensity  (from  Equation  (2.7))  becomes 


h  (u,  v)  =  j  J  d^drii  j  j  d£2dr)2h  (u  -£i,v-  rji)  h*  (u  -  &,v  -  rj2) 

—  OO  —CO  —CO  —CO 

x<4(£i,7i;6,72)  (2.19) 


with  the  mutual  intensity  defined  as 

jg  (6, 7i;  6, 72)  =  < Ug  (6, 71;  t)  u*  (£2,  r]2 ;  t))  .  (2.20) 

The  physical  properties  of  the  two  coherence  extremum  (fully  coherent  and  fully  incoher¬ 
ent)  can  be  exploited  to  define  mutual  intensity.  Considering  coherent  light,  all  the  points 
in  the  field  interfere  with  each  other  (statistically  dependent)  and  it  is  characterized  by  the 
mutual  intensity 

Jg  (6, 71;  6, 72)  =  ug  (6,  Vi)  u;  (6, 72)  (2.21) 
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and  the  resulting  image  intensity  is 


Ii(u,v)  = 


h  (u  -  fi,  v  -  771)  Ug  (fi,  77  1)  diidrjx 
x  h*(u-  £2,  v  -  772)  U*  (6,  V2)  d^2dr)2 


—  OO  —  OO 

CO  CO 


—  CO  —co 

CO  CO 


h  (u  —  £,v  —  rj)  Ug  (f ,  r;)  d^dr] 


—  CO  —CO 

CO  CO 


X 


h{u-£,v-ri)Ug  (i,  r])  d£dr) 


—CO  —CO 
CO  CO 


h  (u  —  i,v  —  r])  Ug  (i,  r])  d^dr] 


—00  —co 


(2.22) 


where  the  relationship  between  the  object  and  image  intensity  is  described  by  the  magnitude 
squared  of  an  amplitude  convolution  relationship  between  the  amplitude  impulse  response 
and  object  optical  field. 

For  incoherent  light,  the  object’s  phasor  amplitudes  are  considered  statistically  inde¬ 
pendent  from  each  other  or,  in  other  words,  the  amplitude  at  one  point  on  the  object  does 
not  affect  the  amplitude  at  a  different  point.  The  mutual  intensity  describing  incoherent 
light  is 

Jg  (fi,  Vi,  6,  V2)  =  kI9  (fi,  771)  5  (£1  -  f2, 771  -  772)  (2.23) 

with  k  being  a  real  constant  and  the  resulting  image  intensity  is 


h  (u,v) 


d£idr]i 


d&dr] 2  x  h  (u  -  £1,  v  -  771)  h*  (u  -  £2,  v  -  772) 


—  OO  — OO 


—  OO  — OO 


X  -6,771  -V2) 


(2.24) 
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and  simplifying  gives 


Ii(u,v)  =  K 


—  OO  —  OO 
CO  CO 


=  K 


h(u  —  £,v  —  rj)  h*  (u  —  £,v  —  rj)  Ig  (f ,  77)  d^dr] 
\h  (u  -  f ,  v  -  77)  |2  Ig  (f ,  7 7)  d£dri 


(2.25) 


—  CO  —co 


where  the  image  intensity  is  a  result  of  an  intensity  convolution  between  the  intensity 
point  spread  function  and  the  object  intensity.  Equation  (2.22)  and  (2.25)  confirm  (2.16) 
and  (2.17)  concerning  the  differences  concerning  intensity  calculations  between  incoherent 
and  coherent  illumination. 

In  LADAR,  it  is  common  to  collect  many  images  to  increase  the  signal  to  noise  ratio 
(SNR)  to  better  enable  detection  and  data  processing.  An  additional  benefit  is  the  partial  co¬ 
herence  of  the  illumination  tends  to  go  from  coherent  to  incoherent  when  averaging  collects 
together.  This  fact  is  due  to  the  many  coherent  images  with  correlated  randomly  varying 
phases  and  amplitudes  combining  to  yield  a  statistically  independent  incoherent  image  map. 
Another  way  to  look  at  why  the  coherent  illumination  goes  to  incoherent  illumination  from 
a  statistical  point  of  view  is  through  the  resulting  probability  mass  function  (PMF)  of  a  par¬ 
tially  coherent  system.  The  PMF  of  a  partially  coherent  system  governing  the  probability 
of  photons  hitting  the  detector  within  one  sampling  interval  is  the  the  negative  binomial 
distribution  given  by  [24] 


Pk(K)  = 


T(K  +  M) 
T(K  +  l)r(M) 


1  M 

-K 

1  K 

1  +  — 

1  +  — 

K 

M 

f>-i  -M 


,  k  =  0,1,  2,  3. 


(2.26) 


where  T  (77)  =  (77  —  1)!  for  any  positive  integer  77,  M  is  the  speckle  parameter  dictating 
the  amount  of  coherence,  and  K  is  the  expected  number  of  photons.  At  the  limits,  M  —  1 
specifies  totally  coherence  and  M  —  00  leads  to  total  incoherence.  In  practice,  all  systems 
fall  somewhere  in  between  the  extremes,  but  assumptions  can  be  made  about  which  end  of 
the  spectrums  dominates.  A  simple  method  to  check  the  coherence  limits  from  statistics  is 
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to  look  at  the  mean  and  variance  of  the  negative  binomial  distribution  given  by 


Unb 


K 


(2.27) 


(2.28) 


As  the  speckle  parameter  M  increases  towards  infinity  at  the  limit,  the  mean  stays  constant, 
but  the  variance  changes  to 


(2.29) 


resulting  in  the  mean  and  variance  being  equal.  This  fact  is  a  characteristic  of  the  Poisson 
distribution  which  has  been  derived  independently  to  characterize  the  probability  of  photon 
hitting  a  detector  given  incoherent  object  illumination  [24].  Figure  2.2  shows  the  effects 
of  speckle  parameter  increase  on  the  negative  binomial  PMF  as  it  approaches  the  Poisson 


PMF. 


The  following  question  still  needs  to  be  answered  explicitly:  Can  the  3D  FLASH 
LADAR  be  considered  to  be  a  result  of  an  incoherent  imaging  process?  Are  many  cubes 
needed  or  just  one?  The  two  analytical  methods  to  provide  convincing  proof  center  around 
obtaining  a  high  speckle  parameter  in  the  partially  coherent  (negative  binomial)  PMF. 

The  first  method  to  attain  a  high  spatial  speckle  parameter  is  to  take  many  indepen¬ 
dent  collects  of  a  particular  remote  scene.  The  speckle  parameter  for  each  collect  is  added 
together  to  yield  a  combined  parameter  which  is  typically  high  enough  to  assume  incoherent 
object  illumination.  For  example,  with  a  mean  number  of  photons  of  50  ( K  =  50),  it  takes 
a  speckle  parameter  of  about  200  for  the  negative  binomial  PMF  to  appear  Poisson.  This 
fact  means  that  even  if  the  light  is  totally  coherent  (M  =  1)  the  resultant  speckle  parameter 
from  summing  the  collects  would  be  sufficient  to  assume  incoherent  object  illumination. 
The  obvious  question  then  arises:  Does  an  operator  have  both  the  time  and  loitering  capa¬ 
bility  to  take  200  collects  of  the  exact  same  scene?  Assuming  the  LADAR  takes  6.7  fi s  to 
take  one  collect  at  a  range  of  1000  meters,  200  collects  would  only  require  1.3  ms  which 
is  a  reasonable  amount  of  time.  The  assumption  that  the  200  collects  consist  of  exactly  the 
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0.06 


K 

(a) 


Figure  2.2:  (a)  This  figure  shows  the  negative  binomial  PMF  for  increasing  values  of  the 

speckle  parameter  at  a  mean  photon  count  of  K  =  50.  As  M  increases,  the 
probability  gets  more  Poisson-like  with  the  main  hump  centered  on  the  mean 
photon  count. 

(b)  This  figure  shows  the  negative  binomial  PMF  at  a  speckle  parameter  of 
M  =  200  and  a  mean  photon  count  of  K  =  50  compared  with  the  Poisson 
PMF  with  the  same  mean. 
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same  scene  is  more  troublesome.  If  an  airborne  platform  is  targeting  in  the  direction  of 
its  velocity,  then  the  consistent  scene  could  be  realized  (if  the  target  isn’t  moving  either). 
However,  as  the  laser  firing  direction  shifts  to  either  side,  the  scenes  are  most  likely  rapidly 
changing  due  to  typical  airborne  platform  speeds.  A  mitigation  to  the  changing  scenes  is 
to  use  available  3D  image  registration  algorithms.  With  respect  to  averaging,  Bayesian 
estimation  attempts  have  been  made  to  mitigate  the  shifts  between  cubes  for  a  particular 
pixel  using  partial  coherent  light  [63].  Since  the  3D  FLASH  LADAR  used  in  this  research 
mounts  on  a  tripod  and  can  easily  obtain  200  or  more  collects  of  a  scene,  the  LADAR  can 
be  used  to  collect  data  with  the  approximation  of  incoherent  object  illumination. 

The  second  method  to  ascertain  if  the  speckle  parameter  is  large  enough  is  from  direct 
calculation.  From  [24],  the  overall  speckle  parameter  M  can  be  defined  as 

M  =  MsMt  (2.30) 


where  Ms  and  M,  are  the  spatial  and  temporal  degrees  of  freedom  respectively.  Given  the 
operating  configuration,  the  area  of  the  detector  Ad  is  smaller  than  the  coherence  area  Ac 
resulting  in  Ms  =  1  [24].  The  area  of  the  detector  is  Ari  =  (100  fi m)2  =  10  nm2  while  the 
coherence  area  Ac  is  defined  by  the  amount  of  coherence  present  in  the  light  given  by  [24] 


Ar  = 


|  n  (Ax,  Ay)  fdAxdAy 


— oo  —  oo 


(2.31) 


with  y,  (Ax,  Ay)  as  the  complex  coherence  factor  that  provides  a  measure  of  the  amount 
of  coherence  between  two  points  and  (Ax,  Ay)  are  the  difference  in  coordinates  between 
two  points  in  the  observation  plane.  In  the  imaging  case,  it  is  shown  that  for  any  incoherent 
source  that 

4,  =  Op-  (2.32) 

A-s 

with  A  as  the  mean  wavelength,  /  as  the  focal  length,  and  As  as  the  area  of  the  incoherent 
light  source.  For  a  circular  aperture,  the  area  of  the  light  source  is  As  =  nr2.  Thus,  the 
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coherence  area  becomes 


„  ((1.55  /im)  (0.3  m))2 

—  n  \2 

7r(l  mm) 


69  nm2 


(2.33) 


and  it  can  be  seen  that  Ms  =  1  due  to  <  Ac  or  10  nm2  <  69  nm2.  The  other  part  of  the 
overall  speckle  parameter  is  the  temporal  degree  of  freedom  Mt  which  is  defined  for  a  light 
beam  with  a  rectangular  power  spectral  density  by  [24] 


A  A 

Mt  =  —  =  j- —  (2.34) 

/Av 

where  A  is  the  pixel  integration  time,  rc  is  the  coherence  time,  and  Av  is  the  bandwidth  of 
the  laser  light.  The  mean  frequency  of  the  laser  light  v  is 


c  3x  108  m/s 

A  1.55  fi m 


194  THz 


(2.35) 


and  assuming  a  bandwidth  of  ±0.05  /mi  gives  a  frequency  bandwidth  Av  =  12.5  THz. 
Considering  an  integration  time  A  =  1  ns,  the  resulting  temporal  degrees  of  freedom  Mt  is 


/\  \  21S 

Mt  =  t - =  1 - - - =  12500.  (2.36) 

lAv  712.5  THz 

Consequently,  the  overall  speckle  parameter  is  M  =  12500  which  is  most  likely  high 
enough  to  assume  incoherent  imaging  by  considering  the  Poisson  distribution  a  valid  ap¬ 
proximation  for  the  negative  binomial  distribution.  This  assumption  would  probably  still 
be  valid  even  if  Av  or  Mt  is  reduced  by  several  orders  of  magnitude. 


2.2  Deconvolution 

With  the  optical  system  able  to  be  represented  by  a  linear  system,  the  attention  turns 
to  the  main  topic  area  in  this  research:  range  estimation  from  object  retrieval  from  data 
observations  of  a  3D  FLASH  LADAR.  The  system  is  modeled  by  a  linear  system  charac¬ 
terized  by  an  impulse  response.  The  observed  data  is  modeled  as  being  generated  from  a 
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convolution  between  the  object  and  impulse  response  corrupted  by  noise.  For  this  research, 
the  object  primarily  consists  of  recorded  amplitudes  and  range  location  of  the  target  under 
interrogation  by  the  3D  FLASH  LADAR.  In  order  to  retrieve  the  object,  the  effects  of  the 
convolution  and  noise  must  be  reversed.  In  other  words,  one  must  deconvolve  the  object 
from  the  impulse  response  while  minimizing  noise  effects.  As  such,  a  review  of  standard 
deconvolution  theory  is  warranted.  The  chosen  model  in  this  research  is  in  units  of  detected 
photons  per  second  while  the  image  intensity  has  only  been  defined  thus  far.  If  ID  denotes 
the  intensity  at  the  detector  (watts/m2),  then  the  following  conversion  results  in  detected 
photons  per  second,  or  mean  photon  flux:  [30] 


(2.37) 


where  A  is  the  cross-sectional  area  of  the  incident  light,  h  is  Planck’s  constant  (6.626  x 
10-34  Joules  •  sec),  and  /  is  the  light’s  frequency.  Substituting  Equation  (2.25)  into  Equa¬ 
tion  (2.37)  gives  the  photons  per  second  at  (it,  v )  in  the  detector  plane  as 


OO  CO 


(2.38) 


where  the  units  would  be  congruent  to  the  mathematical  model  for  the  returned  signal 
presented  in  the  future  sections. 

In  physical  measurements,  noise  mitigation  and  an  unknown  system  impulse  re¬ 
sponse  make  the  problem  more  difficult.  The  system  impulse  response  may  not  be  known 
in  most  operational  ranging  or  imaging  applications.  Thus,  the  process  of  object  retrieval 
is  termed  blind  deconvolution  due  to  the  unknown  system  impulse  response.  In  this  case, 
estimates  of  the  impulse  response  need  to  be  calculated  along  with  the  object  estimates. 

2.2.1  Inverse  Filtering.  If  there  is  no  noise  term  and  the  system  impulse  response 
is  known,  the  deconvolution  can  be  performed  easily  in  the  spatial  frequency  domain.  Note 
that  the  previous  convention  concerning  image  and  object  planes  changes  from  (u,  v)  and 
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(£,  rj)  to  (x,  y )  and  (m,  n)  respectively.  Taking  the  2D  Fourier  transform  of  noiseless  ob¬ 
servations  from  a  3D  FLASH  LADAR,  dk  (x,  y )  in  the  (x,  y )  spatial  dimensions  (ft  is  the 
time  dimension)  results  in 


Dk(fx,  fy)  =  Ok(fx,  fy)H(fx,  fy)  (2.39) 


where  Ok(fx,  fy)  and  H(fx,  fy)  are  the  Fourier  Transform  of  the  object  and  the  system 
impulse  response  respectively.  The  object  can  be  retrieved  by  setting  the  filter,  G,  as  the 
inverse  of  the  Fourier  Transform  of  the  point  spread  function 


ok(m,n)  =  3  1  {Dk(fxJy)G  (fx,  fy)} 
ifDkjfxJy)} 

\  H(fx,fy)  J 

^  ,  \(h(f,.fy)IHf,:fy)\ 

X  H(fXJy)  j 

=  (2.40) 


Conversely,  the  following  highlights  the  severe  limitation  of  inverse  filtering  when  random 
noise  effects  are  introduced: 


Ok(m,  n) 


3 


-1  [DkVxJy) 

X  H(fx ,  fy) 

~l(Ok(fX,fy)H(fx,fy)+Nk(fx,fy 


r-1 


Ok(fx,  fy)  + 


H(fX,fy) 

Nk(fX,fy ) 


H{fx,fy 


(2.41) 

(2.42) 


The  inverse  filter  solution  will  be  skewed  to  the  degree  that  the  impulse  response  amplifies 
the  noise  term  which  can  be  significant.  This  noise  amplification  is  a  primary  driver  towards 
other  solutions  based  on  minimizing  the  effects  of  noise. 
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According  to  [37],  the  Wiener  filter  minimizes  the  mean  squared  error  between  the 
real  object  and  the  estimated  object,  E[(o  —  o)2],  resulting  in  the  following  functional  form 


G(fX,  fy ) 


(2.43) 


I  H(fx,  fy)\2  +  Snn{fx,  fy)/Su(fx ,  fy) 


where  *  is  the  conjugate  operator  and  Snn  and  Stt  are  the  power  spectra  of  the  noise  and 
signal  respectively.  The  resulting  estimate  for  the  object  is 


(2.44) 


Examining  this  final  form  is  enlightening  to  how  the  filter  handles  certain  noise  situations. 
When  the  noise  spectrum  is  zero  or  dominated  by  signal,  the  filter  simplifies  to  the  inverse 
filter.  When  the  noise  power  is  severe  or  the  signal  level  is  low  at  some  frequencies,  the 
filter  approaches  zero  attenuating  these  frequencies  with  high  noise  power. 

2.2.2  Iterative  Algorithms.  Iterative  deconvolution  techniques  also  exist  to  in¬ 
clude  the  Richardson-Lucy  and  error  minimization  algorithms  which  are  useful  when  data 
models  are  complex  or  non-linear.  From  [62],  the  Richardson-Lucy  algorithm  was  de¬ 
veloped  to  be  an  approximate  deconvolution  to  recover  the  object  W  from  the  degraded 
noiseless  image  H  =  W  0  S  with  %  as  the  convolution  operator  and  S  as  the  point  spread 
function.  Temporarily  adopting  notation  from  [62],  the  problem  is  constructed  based  on 
Bayes  theorem  given  by  [60] 


(2.45) 
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where  P(W\H)  is  the  conditional  probability  of  W  given  H  (also  called  the  a  posteriori 
density),  P(H \W)  is  the  conditional  probability  of  H  given  W,  P(W)  is  the  marginal  prob¬ 
ability  of  W  (also  called  the  a  priori  density),  and  P(H )  is  the  marginal  probability  of  H. 

The  subscripts  j  and  q  correspond  to  pixel  locations  with  W  =  W3  and  H  =  Hq 

j  i 

equalling  the  value  for  the  entire  object  and  degraded  image  arrays  respectively.  The  prior 
probability  can  be  defined  by  [60] 


(2.46) 


and  by  combining  Equations  (2.45)  and  (2.46)  results  in  the  following  equation  [62] 


(2.47) 


Noting  that  the  desired  solution,  P  (Wj),  is  also  on  the  right-hand- side  of  the  equation  and 
is  not  a  function  of  the  summation,  a  common  practice  is  to  make  an  initial  guess  and  set 
up  the  iterative  updates  as 


(2.48) 


Reduction  of  Equation  (2.48)  is  still  necessary  due  to  being  in  terms  of  probability.  This 
equation  is  changed  so  that  it  uses  actual  variable  values  rather  than  probability.  Using  the 
laws  of  probability  and  the  conservation  of  energy,  the  probabilities  can  be  reformed  into 


P(Wj)  =  Wj/W, 

P  (Hq)  =  Hq/H  =  Hq/W , 


and  P  (Hq\Wj)  =  P(Sm)  =  S3JS. 


(2.49) 
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Consequently,  Equation  (2.48)  can  be  reduced  to 


Wj>r+1  =  W, 


v-  SjrEz 


hr 


2. 


E  Si,zWi, 


(2.50) 


which  represents  the  final  form  of  the  Richardson-Lucy  object  recovery  algorithm.  One 
weakness  of  this  algorithm  is  its  lack  of  proven  convergence.  In  practice,  however,  the 
iterations  provide  the  perfect  solution  in  the  noiseless  case  and  an  improved  solution  with 
noisy  data. 

From  [8],  the  last  reviewed  method  of  deconvolution  involves  using  a  cost  function 
and  minimizing  it  with  respect  to  the  data  and  the  true  image.  The  cost  function  is  defined 
by 

M  N 

C  =  ^2^(d(x,y)  -  i(x,y))2  (2.51) 

x=l  y= 1 

where  the  data  equals  d  (x,  y)  —  i  (x,  y)+n  (x,  y)  where  n  (x,  y)  is  the  signal  independent, 
additive  noise  and  the  true  image  is  defined  as 

M  N 

i  (x,  y)  =  o  (m,  n)  h  (x  —  m,  y  —  n)  (2.52) 

?7i=l  n= 1 


with  o  (m,  n )  as  the  object  and  h  (m,  n )  as  the  point  spread  function.  In  order  to  minimize 
the  cost  with  respect  to  the  unknown,  the  derivative  of  the  cost  function  is  taken  with  respect 
to  the  object  with  the  result  set  to  zero.  The  solution  is  obtained  by  solving  this  equation  for 
the  object.  Thus,  the  derivative  of  Equation  (2.51)  is  taken  with  respect  to  a  single  object 
pixel  (m0,  n0)  and  set  to  zero 


dC 

do  (m0,  na) 


M  N 

-2  (d  y) 

x=l  y=  1 


(%,y)) 


di  (x,  y) 
do  (m0,  nQ) 


0. 


(2.53) 
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The  partial  derivative  of  the  image  is 


di  (x,  y) 
do  ( m0,n0 ) 


- - -  [o  (m0,  n0)  h(x  -  m0,y  -  n0)\ 

do  (m0,  n0) 

h(x  -  m0,y  -  nQ) 


(2.54) 


giving  the  resulting  expression  and  reduction 

M  N 

EE  (d  (x,y)  —  i  (x,y))h  (x  —  m0,y  —  n0)  =  0 

x=l  y=  1 

M  N 

E  E  d  (,x>  y)h  ( x  -  mo,  y  -  na) 

x=l  y=  1 _  _ 

~~M  N  _ 

E  E  *  (x,  y)h  ( x  —  m0,  y  —  na) 

x=l  y=  1 


(2.55) 


Using  reasoning  similar  to  [62],  the  object  is  then  multiplied  on  both  sides  of  Equa¬ 
tion  (2.55)  giving  the  final  form  of  the  object  recovery  as 


M  N 

E  E  d  (x,  y)h  {x  —  m0,  y  —  na) 

onew  (m0,  n0)  =  oold  (m0,  na)  -  (2.56) 

E  E  iold  (xi  y)h  ( x  -m0,y-  nQ) 

x=l  y= 1 


M  N 

with  iold  (x,  y)  =  E  E  °°ld  (m>  n)  h  (x  —  m,  y  —  n) .  An  acceptable  stopping  point  can 

m= 1 n=l 

be  (1)  minimal  change  from  the  previous  iteration  or  (2)  the  appropriate  amount  of  image 
noise  in  the  estimated  image  based  on  prior  knowledge  of  the  noise  source. 


2.3  Maximum  Likelihood 

Maximum  Likelihood  (ML)  estimation  can  be  used  for  a  data  model  that  includes  no 
blurring  function  because  the  model  then  implicitly  assumes  no  coupling  between  pixels. 
The  ML  method  can  then  operate  on  one  pixel  at  a  time.  Another  means  of  estimating 
parameters  (e.g.  object,  received  amplitude,  range  to  target)  is  to  employ  ML  estimation 
using  the  observation  statistics  to  form  a  likelihood  expression.  Lrom  [84],  the  ML  solution 
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is  the  outcome  of  a  ML  analysis  where  the  estimate,  a,  maximizes  the  likelihood  function, 
L(A),  or 

a  —  argrnax  L{A)  (2.57) 

where  the  parameter  A  can  either  be  a  single  or  vector  variable.  Considerations  of  max¬ 
imum  likelihood  estimation  include  the  uncertainty  that  a  unique  ML  solution  exists  and 
local  maximums  in  the  likelihood  function. 

One  way  to  view  the  ML  solution  is  as  a  special  case  of  Maximum  a  Posterior  (MAP) 
estimation  with  the  prior  distribution  being  a  uniform  distribution.  MAP  estimation  is 
Bayesian  based  and  starts  with  Bayes  Theorem.  Recall  that  Bayes  Theorem  relates  the 
conditional  and  marginal  probabilities  of  events  A  and  B  with  B  having  a  non-zero  proba¬ 
bility.  The  equation  for  Bayes  Theorem  is  defined  again  as  [60] 


Pa\b(A\B) 


Pbla(B\A)Pa(A ) 

Pb(B) 


(2.58) 


where  P(A\B)  is  the  conditional  probability  of  A  given  B  (also  called  the  a  posteriori  den¬ 
sity),  P(B\A)  is  the  conditional  probability  of  B  given  A,  P(A)  is  the  marginal  probability 
of  A  (also  called  the  a  priori  density),  and  P(B)  is  the  marginal  probability  of  B.  Bayes 
theorem  calculates  the  probability  of  event  A  occurring  given  observing  B.  Maximizing 
Equation  (2.58)  is  mathematically  equivalent  to  maximizing  the  natural  log  resulting  in 


ln[pa\b(A\B)]  =  In pb\a(B\ A)  +  lnpa(A)  -  In pB{B).  (2.59) 


The  MAP  estimate  is  found  by  taking  the  derivative  of  Equation  (2.59),  setting  it  equal  to 
zero,  solving  for  A  given  by 


<91n(p(A|T?))  <91n(p(f?|A))  <91n(p(A))  <91n(p(f?)) 

BA  =  dA  +  dA  dA 


(2.60) 
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where  aill(.^ 1 1  =  0  due  to  no  dependence  on  A.  The  final  form  of  the  MAP  estimator  is 
then 


Carnap  arg  max  (lnp6|a  (B\A)  +  In pa  ( A ))  .  (2.61) 

When  the  prior  probability  pa(A)  is  unavailable  or  not  postulated,  it  can  be  assumed  that 
the  prior  probability  can  be  described  as  a  uniform  RV.  Thus,  pa(A)  has  no  dependence  on 
A  either  and 


<91n(p(A)) 

(M 


0 


(2.62) 


resulting  in  the  ML  solution  of 


amap  =  arg  max  (hip  (B\Aj) .  (2.63) 

A  maximum  likelihood  technique  is  used  for  single  pixel  range  estimation  in  Section  4.2. 


2.4  Generalized  Expectation  Maximization 

Traditional  linear  maximum  likelihood  efforts  do  not  suffice  to  estimate  target  range 
given  the  unknowns  (amplitude,  target  range,  PSF,  and  pixel  bias)  in  the  statistical  model 
from  Equation  (2.69).  More  powerful  object  estimation  methods  like  the  Generalized  Ex¬ 
pectation  Maximization  (GEM)  algorithm  must  be  employed  due  to  the  coupled  unknowns 
which  will  be  covered  in  the  next  section.  While  the  final  goal  is  to  estimate  range,  a  dif¬ 
ferent  tactic  is  employed  due  to  the  difficulty  in  having  the  target  range  term  residing  in 
the  exponential.  Consequently,  the  unknowns  in  the  estimation  process  are  the  target  am¬ 
plitude,  target  pulse  shape  (or  object),  and  PSF.  With  the  pulse  shape  now  as  an  unknown, 
it  is  much  simpler  to  use  the  GEM  to  find  maximum  likelihood  solutions.  Once  the  max¬ 
imum  likelihood  solution  for  the  object  or  pulse  shape  is  found,  a  correlation  operation 
between  the  estimated  pulse  shape  and  a  reference  pulse  shape  determines  the  estimated 
target  range.  A  full  description  of  the  algorithm  will  be  given  in  the  subsequent  paragraphs. 
First,  the  GEM  solutions  for  the  unknown  parameters  must  be  found.  However,  a  closed 
form  solution  for  the  EM  algorithm’s  maximization  step  is  intractable.  Consequently,  the 
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GEM  algorithm  goal  is  to  modify  the  EM  structure  such  that  the  likelihood  is  incremen¬ 
tally  increased  rather  than  globally  maximized  as  in  the  EM  algorithm.  This  incremental 
increase  in  the  likelihood  simplifies  the  maximization  step  allowing  unknown,  non-random 
parameter  estimation. 

In  the  case  of  blind  deconvolution,  the  EM  algorithm  can  be  implemented  to  esti¬ 
mate  the  object,  point  spread  function,  range,  and/or  amplitude.  This  algorithm  is  a  another 
method  to  perform  maximum  likelihood  estimation  whereby  the  solution  is  found  by  us¬ 
ing  unobserved  data  (complete)  rather  than  the  observed  data  (incomplete).  Although,  the 
maximum  likelihood  solution  is  not  always  guaranteed  as  a  result  from  the  EM  algorithm. 
Pertaining  to  the  unobserved  data,  it  may  be  necessary  because  the  regular  maximum  like¬ 
lihood  solution  may  be  analytically  prohibitive.  The  EM  algorithm  uses  the  reduced  com¬ 
plexity  of  the  complete  data  problem  to  perform  maximum  likelihood  estimation. 

According  to  [54],  the  EM  algorithm  is  composed  of  two  steps.  The  first  step  (E- 
Step)  is  to  find  Q :  the  expected  value  of  the  desired  variable  given  the  latest  parameter 
values  or 

Q  (*;  *(fc))  =  £*(*)  {In  Lcd  (1®0  |  y}  (2.64) 

where  is  the  vector  of  unknown  parameters,  k  is  the  iteration,  LCd{^)  is  the  complete 
data  likelihood,  and  the  expectation  conditioned  on  the  incomplete  data  y.  Complete  data 
can  be  viewed  as  the  unobserved  variables  (fabricated  or  not)  used  to  simplify  the  problem. 
Incomplete  data  is  usually  the  observed  data.  The  second  step  (M-Step)  is  to  maximize  this 
expected  value  with  respect  to  the  unknown  parameters,  \E,  by  choosing  to  maximize 
Q  (’E;  ^(fc))  or 

Q  (^(fe+1);  >  Q  (¥;  ¥(fc) )  .  (2.65) 

for  all  unknown  parameters  in  \P.  The  EM  algorithm  is  advantageous  due  to  the  guarantee 
of  increasing  the  likelihood  with  each  iteration  and,  in  most  cases,  eventually  converging  on 
the  maximum  likelihood  solution.  As  proven  by  [16],  the  incomplete-data  log-likelihood 
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function  increases  with  each  iteration 


L  (^(fc+1))  >  L  (^(fc)) 


(2.66) 


and  the  EM  algorithm  converges  to  local  or  global  maximum. 

As  is  the  case  in  this  research  where  the  maximization  over  all  unknown  parameters 
is  difficult  or  doesn’t  exist  in  a  closed  form,  an  incremental  EM  algorithm  is  used  called 
the  generalized  expectation  maximization  (GEM)  where  the  goal  is  to  simply  increase  the 
likelihood  at  each  iteration  without  finding  the  maximum  parameter  value.  A  GEM  requires 
that  the  likelihood  be  improved  and  not  maximized  such  that 


(2.67) 


If  Equation  (2.67)  holds  for  every  iteration,  it  has  been  shown  that  the  likelihood  is  in¬ 
creased  with  every  iteration  or  [54] 


L  (^(fc+1))  >  L  (^(fc)) 


(2.68) 


and,  if  bounded,  the  GEM  sequence  converges  to  a  local  maximum  due  to  the  monotonicity 
of  the  algorithm.  The  GEM  algorithm  will  be  implemented  on  simulated  and  experimental 
data  in  Chapter  V  to  show  that  object  recovery  improves  range  estimation. 

2.5  3D  FLASH  LADAR  Data  Model 

This  section  describes  the  physical  3D  FLASH  LADAR  model.  To  increase  read¬ 
ability,  the  model  is  defined  in  this  chapter  due  to  parts  of  Chapter  IV  and  all  of  Chapter  V 
using  this  particular  model.  (Other  sections  in  Chapter  IV  and  all  of  Chapter  VI  use  a 
different,  simplified  observation  model  to  allow  for  relatively  uncomplicated  mathemati¬ 
cal  expressions  and  for  concept  investigation.  The  changes  in  model  definition  are  clearly 
identified.) 
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Figure  2.3:  (a)  3D  view  of  LADAR  system  model  in  Cartesian  coordinates  with  each 

data  cube  having  dimensions  of  30  x  30  x  20  corresponding  to  pixel  x  pixel 
x  time  sample.  The  variable  d(tk)  corresponds  to  the  kth  receiver  detected 
range  slice  image  with  k  e  [1, ...,  N]  and  N  =  20. 

(b)  Another  view  of  the  3D  FLASH  LADAR  operation.  The  N  number 
of  samples  are  meant  to  depict  the  available  target  information  that  the  2D 
range  images  (slices)  would  collect.  The  assumed  time  separation  between 
the  range  images  is  2  nanoseconds  closely  corresponding  to  the  3D  LADAR 
system  used  for  experimental  collects. 
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Simple  Model 


d(x,y) 


High  Fidelity  Model 


Figure  2.4:  For  a  given  range  slice,  this  diagram  shows  the  propagation  of  the  object 
through  optical  system  to  the  observation.  Definitions:  o  is  the  object,  h  is 
the  PSF,  B  is  the  pixel  bias,  n  is  the  noise,  and  d  is  the  observation.  The 
simple  signal  model  is  used  in  previous  3D  FLASH  LADAR  research  such 
as  [9],  [39],  and  [55].  The  high  fidelity  model  is  used  in  Chapters  V  and  VI. 


Figure  2.3  shows  the  sensor  operation  resulting  in  a  data  cube  of  spatial  and  range 
information.  In  simple  terms,  the  LADAR  laser  transmits  a  pulse  and  the  LADAR  detector 
array  receives  an  attenuated,  time-delayed  version  of  the  transmitted  pulse.  Each  detec¬ 
tor  receives  a  version  of  the  waveform  shape  sampled  according  to  the  range  gate.  Thus, 
models  can  take  advantage  of  this  fact  and  perform  range  estimation  on  a  per  pixel  basis. 
Referring  to  Figure  2.4,  previous  research  has  assumed  the  simple  model  [9],  [39]  where 
the  spatial  impulse  response  was  a  Dirac  delta  function.  This  definition  meant  there  were 
no  interactions  between  adjacent  pixels.  However,  the  research  in  this  dissertation  adopts 
the  high  fidelity  model  since  it  is  more  accurate  concerning  pixel  spatial  interactions.  The 
limitation  of  the  simplistic  model  and  adaptation  of  the  higher  fidelity  model  is  the  catalyst 
of  the  material  in  Chapters  V  and  VI. 

In  order  to  simplify  the  geometry  and  the  mathematics,  assumptions  are  made  about 
the  model  to  include:  (i)  target  is  perpendicular  to  the  transmitter,  (ii)  target  is  in  the  far-held 
of  the  receiver,  (iii)  target  is  Lambertian,  (iv)  circular  optics  are  in-focus,  (v)  monostatic 
RADAR  operation,  (vi)  the  waveform  is  a  symmetric  Gaussian  pulse,  (vii)  each  pixel  from 
the  detector  array  has  an  individual  waveform  associated  with  it,  and  (viii)  the  range  slices 
of  the  data  cube  are  registered.  Other  pulse  shape  models  are  available  include  an  asym¬ 
metric  Gaussian,  a  truncated  negative  parabolic,  or  some  hybrid  of  a  Gaussian  and  negative 
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parabolic.  A  symmetric  Gaussian  is  chosen  for  notation  purposes,  but  an  asymmetric  ver¬ 
sion  is  easily  defined  with  different  pulse-widths  for  the  leading  and  trailing  edges. 

Considering  a  3D  FLASH  LADAR  sensor  with  statistically  independent  samples 
dominated  by  shot-noise  [9],  [39],  the  PMF  of  the  observed  photons,  djk(x,y),  incor¬ 
porating  all  cubes  ( j  G  [1, ...,  J]),  range  samples  (k  G  [1, ....  A']),  and  detector  pixels 
(x  G  [1,  X\  ,  y  G  [l,...,y])is 


P  [ Djk  (x,  y)  =  djk  (x,  y) ;  Vj,  k,  x,  y]  = 

TT  [ijk  {x,  y)  +  B  (x,  y)]d^x'v)  exp  {-  [ijk  (x,  y)  +  B  (x,  y)}} 
ill,  dlk(x,y)i 

where  the  mean  signal  is  ijk  (x,  y)  +  B  (x.  y)  where  B  (x,  y)  is  the  pixel  bias  and  the  blurry, 
non-noisy  signal  ijk  (x,  y)  is  defined  by 

M  N 

ijk  (x,  y)  =  ^2  ^2  °k  ("A n)  hj  (x  -  m,y  -  n)  (2.70) 

771=1  71=1 

where  the  object  ok  (m,  n )  is  defined  at  the  object  plane  with  coordinates  (m  G  [1, ...,  M] 
and  n  G  [1, ....  A']).  The  (x,  y)  and  k  variables  correspond  to  a  pixel  in  the  detector  array 
and  to  the  returned  signal  time  of  arrival  respectively.  The  time  of  arrival  is  computed  based 
on  the  time  from  laser  pulse  transmission  to  photon  detection.  This  assumption  may  require 
cube  registration  due  to  the  possibility  of  moving  targets,  moving  sensor  platform,  or  inter¬ 
cube  timing  errors.  Incorporating  contributions  from  light  propagation,  optical  abberations, 
and  atmospheric  blurring,  the  intensity  point  spread  function  (PSF)  hj  (x,  y)  is  constant 
within  a  single  cube  while  changing  across  cubes.  In  this  research,  the  PSF  is  considered 
constant  within  a  single  cube  since  collection  times  spans  under  forty  nanoseconds  and 
any  time-dependent  effects  would  be  essentially  frozen.  In  addition,  the  pixel  bias  B  (x,  y) 
is  constant  between  cubes  as  well  as  within  a  single  cube  due  to  the  pixel’s  unchanging 
physical  material  and  response  to  incident  light  (ambient  radiation  is  assumed  negligible). 

Every  pixel  in  the  detector  array  records  a  time-delayed  and  attenuated  version  of  the 
transmitted  pulse.  The  physical  outgoing  pulse  shape  of  a  3D  FLASH  LADAR  is  either 
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Gaussian,  negative  parabolic,  or  some  hybrid  of  the  two.  The  object  can  be  defined  by  an 
amplitude  term  and  a  pulse  shape  or 

Ok  (m,n)  =  A(m,n)  pk  (m,n) .  (2.71) 


Assuming  a  Gaussian  transmitted  pulse,  the  object  is 


Ok  (m,  n) 


A  (m,  n ) 


-  (4 


2 R  (m,  n)/cf 


(2.72) 


where  A  (m,  n)  is  the  object  amplitude,  aw  is  the  waveform  standard  deviation,  tk  is  the 
time  variable,  c  is  the  speed  of  light,  and  R  (m,  n )  is  the  range  to  the  target.  If  a  negative 
parabolic  waveform  model  is  desired,  the  object  is  defined  by 


ok  ( m ,  n)  =  A  (m,  n ) 


(2 R  (m,  n)  —  tkc )' 
( cpw )2 


rect 


2 R  (m,  n)  —  tkc 
2  cpw 


(2.73) 


where  2 pw  is  the  pulse  width  and  rect  is  the  rectangle  function  defined  by 


0,  if  |x|  >1/2 


rect  (x) 


<1/2,  if  |x|  =  1/2 


1,  if  |x|  <  1/2. 


(2.74) 


Although,  for  simplicity  and  ease  of  differentiation,  this  research  adopts  the  Gaussian 
model.  For  military  targeting  or  navigation,  range  to  target  (located  in  the  object  term) 
is  the  desired  unknown  variable.  In  attempting  to  perform  range  estimation,  a  range  term 
is  not  explicitly  in  the  model,  but  it  is  buried  within  the  object,  ok  (m,  n),  term  given  by 
Equation  (2.72)  or  (2.73).  If  the  object  were  exactly  known,  the  target  range  could  be  then 
extracted  from  the  object  by  peak  detection  methods.  This  statement  presents  the  ideal  sit¬ 
uation  that  Chapter  V  attempts  to  create  with  an  object  degraded  by  spatial  blurring  and 
noise  sources. 
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Given  the  LADAR’s  3D  nature,  it  is  important  to  discern  the  formation  of  range  slice 
images  shown  in  Figure  2.3  versus  the  pixel  waveform  definition  from  (2.72).  The  range 
slice  images  are  formed  at  a  particular  time  by  a  spatial  convolution  between  the  original 
scene  and  the  system’s  impulse  response.  The  original  scene’s  amplitude  is  variable  at  each 
time  instant  due  to  target  roughness  and  Gaussian  shaped  transmitted  pulse.  Therefore, 
the  returned  amplitude  changes  for  each  range  image  formation  operation.  Considering 
atmospheric  turbulence,  the  system’s  impulse  response  is  assumed  constant  for  each  [1,  N] 
range  image  due  to  the  short  duration  of  the  data  cube  collection  (forty  nanoseconds)  [24]. 
Conversely,  the  pixel  waveform  definitions  from  Equation  (2.72)  define  each  pixel’s  un¬ 
blurred  and  non-noisy  received  signal  where  the  model  assumes  only  one  target  per  pixel. 
The  range  estimation  process  estimates  the  target’s  range  for  every  pixel.  The  following  is 
a  concise  explanation  of  the  difference  between  data  generation  and  range  estimation:  the 
simulation  forms  3D  LADAR  data  cubes  in  the  spatial  domain  while  the  range  estimator 
operates  in  the  range  (time)  domain.  Also,  as  will  be  discussed  later,  image  deblurring 
operates  spatially  like  the  image  formation  process. 

Following  [25],  a  transfer  function  describes  the  LADAR’s  effect  on  the  target  return 
assuming  the  system  is  linear  and  spatially  invariant.  The  transfer  function  in  optics  is  called 
an  optical  transfer  function  (OTF).  If  only  considering  the  effects  of  the  optical  components, 
the  OTF  is  diffraction-limited  because  the  only  way  to  increase  performance  would  be  to 
build  better  optical  components.  Otherwise,  optical  diffraction  theory  bounds  the  system 
performance. 

While  not  the  main  focus  of  this  research,  it  is  important  to  understand  that  3D 
FLASH  LADAR  operational  use  may  encounter  periods  of  considerable  atmospheric  tur¬ 
bulence  that  would  modify  the  system  OTF.  As  long  as  the  imaging  scenario  stays  within 
the  isoplanatic  angle,  the  PSF  can  still  be  considered  spatially  invariant  which  is  a  prereq¬ 
uisite  to  this  deblurring  technique  [51].  Given  this  condition  holds  true,  the  OTF  is  then  a 
function  of  the  diffraction-limited  OTF  and  the  average  OTF  resulting  from  the  atmosphere. 
Considering  a  substantial  target  distance,  there  could  be  atmospheric  distortion  and  would 
manifest  itself  by  modifying  the  diffraction-limited  OTF  to  form  an  overall  OTF  [24].  Ne- 
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Figure  2.5:  An  example  of  a  diffraction  limited  OTF.  This  OTF  was  generated  using  the 
parameters  from  this  research. 

glecting  pixel  integration  effects,  the  form  of  the  overall  OTF,  H(fx,  /,/),  could  be 


H{fX,fy)=H0{fx,fy)HA{fx,fy) 


(2.75) 


where  (fx,  fy)  are  spatial  frequency  variables,  H0(fx,  fy)  is  the  diffraction-limited  OTF,  and 
HA(fx,  fy)  is  the  short-exposure  average  OTF  due  to  atmospheric  turbulence.  The  form  of 
Ha  is  [24]  [67] 


with  v  =  sj fx  +  fy,  A  the  mean  wavelength,  /  the  optic’s  focal  length,  rQ  as  the  atmo¬ 
spheric  coherence  diameter  or  Fried’s  seeing  parameter,  and  D  is  the  aperture  diameter. 
With  the  OTF  defined  as  the  inverse  Fourier  Transform  of  the  PSF,  Figure  2.5  shows  a  two- 
dimensional  representation  of  the  simulation’s  diffraction-limited  OTF.  Using  centered  ID 
cutouts,  Figure  2.6  shows  the  effect  of  the  atmosphere  on  the  OTF  whereby  the  atmosphere 
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amplitude 


Figure  2.6:  Cut-outs,  (0,  fy),  of  different  OTFs  to  include  diffraction-limited,  atmo 
sphere,  and  overall.  The  degradation  in  the  overall  OTF  caused  by  the  at 
mosphere  is  evident  with  higher  spatial  frequencies  blocked. 
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degrades  the  overall  OTF  by  narrowing  the  amount  of  spatial  frequencies  the  system  can 
pass.  This  truncation  of  spatial  frequencies  causes  high  frequency  details  in  the  range  slice 
image  (i.e.  sharp  comers,  fine  lines,  etc.)  to  be  lost.  The  narrowing  of  the  OTF  in  the  spa¬ 
tial  frequency  domain  leads  to  a  widening  of  the  PSF  in  the  spatial  domain.  This  widening 
causes  increased  pixel  mixing  due  the  the  convolution  nature  of  the  system.  The  result¬ 
ing  received  waveform  is  further  deviated  from  the  idealized  received  waveform  in  (2.72). 
Blind  deconvolution  methods  in  Chapter  V  would  effectively  estimate  any  additional  at¬ 
mospheric  blurring  as  long  as  the  mode  of  operation  remained  in  the  isoplanatic  angle  (i.e. 
spatially  invariant)  [25]. 
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Figure  2.7 :  Organizational  chart  for  the  literature  review.  The  review  is  broken  down 

into  the  following  sections:  3D  FLASH  LADAR  Hardware  and  Applica¬ 
tions  -  Section  2.6.1.  3D  FLASH  LADAR  Post-Processing  -  Section  2.6.2. 
Blind  Deconvolution  -  Section  2.6.3.  CRB  and  Parameter  Optimization  - 
Section  2.6.4. 

2.6  Previous  Research 

This  section  contains  the  literature  review  of  publications  relating  to  hardware  devel¬ 
opment  and  post-processing  of  3D  FLASH  LADAR  data.  The  background  review  provides 
a  treatment  of  several  important  topics:  3D  FLASH  LADAR  hardware  development  and 
applications,  3D  LADAR  post-processing  algorithms,  LADAR  range  estimation,  general 
blind  deconvolution  theory  and  applications,  lucky  imaging,  and  3-D  image  registration. 
Seminal  papers  are  reviewed  as  well  as  appropriate  recent  publications.  For  easy  reference, 
Figure  2.7  shows  the  literature  review  organization. 

LADAR  theoretical  development  in  the  past  10-20  years  has  concentrated  on  3D 
scanning  LADAR  systems  almost  exclusively  because  it  was  the  only  3D  LADAR  avail¬ 
able.  3D  FLASH  LADARs  are  a  relatively  new  development,  explaining  the  lack  of  publi¬ 
cations  compared  to  more  mature  technologies.  The  current  3D  FLASH  LADAR  literature 
spans  from  hardware  development  to  applications  to  post-processing.  The  post-processing 
papers  consider  important  algorithms  enabling  improved  range  estimation,  feature  extrac¬ 
tion,  foliage  penetration,  world  modeling,  mapping,  and  navigational  aiding. 

The  papers  that  do  take  advantage  of  the  unique  properties  of  a  3D  FLASH  system 
use  a  simplified  data  model.  The  spatial  convolution  nature  between  the  object  plane  in¬ 
tensity  and  the  detector  plane  intensity  is  not  accounted  for,  leading  to  errors  in  parameter 
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estimation.  There  is  a  gap  in  the  literature  considering  the  spatial  effects  of  a  3D  LADAR 
system  because  the  scanning  systems  simply  don’t  see  the  effects  of  the  spatial  impulse 
response.  Scanning  LADARs  don’t  operate  fast  enough  spatially  or  have  a  wide  enough 
field-of-view  (FOV)  to  observe  the  blurring  effects  of  the  spatial  impulse  response  like  the 
FLASH  systems  do. 

2. 6. 1  3D  FLASH  LADAR. 

2.6. 1.1  Hardware  Developments.  Although  new  advances  LADAR  hard¬ 
ware  development  is  not  the  focus  of  this  research,  it  is  prudent  to  know  about  not  only  the 
hardware  used  in  this  research,  but  also  other  state-of-the-art  3D  LADARs.  Understand¬ 
ing  where  the  technology  stands  and  some  of  the  details  will  lend  an  appreciation  for  the 
uniqueness  and  potential  of  the  3D  FLASH  LADARs.  The  advances  made  in  the  LADAR 
hardware  have  increased  capability,  but  have  created  additional  issues  that  need  mitigation. 

Based  on  work  from  [64],  [65],  and  [66],  the  enabling  technology  allowing  3D 
FLASH  LADAR  to  be  realized  culminated  in  2004  with  the  development  of  a  focal  plane 
array  (FPA)  capable  of  collecting  a  series  of  two  dimensional  (2D)  images  of  a  remote 
scene  at  varying  depths  from  a  single  laser  pulse  [81].  The  modelling  in  this  research  is 
based  on  this  hardware.  Additionally,  this  particular  3D  FLASH  LADAR  system  will  be 
used  for  experimental  data  collection  in  the  future.  The  novel  hardware  design  using  de¬ 
tector  material  made  of  either  InGaAs  PIN  or  Avalanche  Photodiodes  (APD),  along  with 
the  data  acquisition  board  called  the  Readout  Integrated  Circuit  (ROIC),  allows  for  rapid 
data  collection  in  the  range  dimension  with  each  pixel  able  to  digitally  sample  the  returned 
waveform.  The  ROIC  permits  this  rapid  range  sampling  with  a  bank  of  capacitors  behind 
each  pixel  capable  of  operating  on  the  nanosecond  scale.  Of  note,  a  similar  LADAR  de¬ 
veloped  with  the  same  goals  is  summarized  in  [29].  The  only  noticeable  difference  in  this 
LADAR  was  that  it  uses  HgCdTe  APD  detector  technology  exclusively. 

The  FLASH  LADAR  is  considered  an  improvement  over  scanning  LADARs  when 
considering  all  the  scene’s  information  is  collected  in  one  shot  and  that  there  is  no  need  for 
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pixel  registration.  Additionally,  the  3D  FLASH  LADAR  is  “eye-safe”  because  it  operates 
in  the  short  wave  infrared  (SWIR)  regime  (beyond  1.4  microns).  It  has  been  shown  that  the 
selected  detector  materials  perform  well  at  this  wavelength  with  respect  to  both  quantum 
efficiency  and  electrical  bandwidth.  Also,  there  are  substantial  cost  and  weight  savings 
given  that  a  mechanical  steering  mechanism  is  not  needed  like  in  the  scanning  systems. 

While  there  are  obvious  benefits,  there  are  several  drawbacks  to  the  system  that  need 
addressing  in  future  hardware  upgrades.  There  are  a  limited  number  of  range  samples 
available  for  each  transmitted  signal.  Essentially,  there  is  a  limit  to  the  time  the  “shutter” 
can  be  open.  In  one  operating  mode,  this  limits  the  operator  to  know  where  the  target  is  a 
priori  to  within  several  meters.  This  limitation  is  not  an  issue  in  the  laboratory,  but  will  need 
to  be  addressed  for  operational  use  either  in  hardware  upgrades  that  solve  the  problem  or 
by  CONOPS  (Concept  of  Operations).  For  example,  another  sensor  could  roughly  locate 
the  target  and  pass  that  rough  location  to  the  3D  FLASH  LADAR  to  fine  tune  the  range 
measurements.  Another  issue  mitigated  in  [73]  is  pixel  coupling  occurring  throughout  the 
detector  array  caused  by  a  time-dependent  gain  variation.  Finally,  as  mentioned  before, 
spatial  impulse  response  effects  are  now  evident  in  the  data  and  are  the  primary  focus  of 
this  prospectus. 

Advances  in  technology  like  the  AFRL  3D  FLASH  LADAR  are  an  example  of  hard¬ 
ware  improvements  opening  up  fields  of  research  not  otherwise  considered.  Evolving  tech¬ 
nology  from  scanning  to  FLASH  LADARs  will  vastly  increase  operational  capabilities  and 
pave  the  way  for  future  applications. 

Other  efforts  to  produce  3D  FLASH  LADAR  hardware  have  succeeded  as  well. 
In  [31],  advances  in  detector,  electric  circuitry,  and  laser  transmitter  technology  are  dis¬ 
cussed  with  the  capability  to  capture  an  entire  3D  scene  in  one  transmitted  pulse.  The  ad¬ 
vances  are  similar  to  [81]  with  some  minor  differences:  (1)  using  the  APD  in  Geiger-mode 
due  to  laser  compatibility  and  size  and  power  requirements  and  (2)  the  capture  circuitry 
is  CMOS-based  resulting  in  a  0.5  nanosecond  timing  resolution.  This  timing  resolution 
corresponds  to  range  information  (i.e.  taking  a  picture)  every  15  cm  (30  cm  for  the  AFRL 
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3D  FLASH  LADAR).  A  key  point  in  the  paper  is  the  huge  benefit  of  employing  a  photon¬ 
counting  3D  LADAR  with  APD  detectors  versus  a  CCD  camera  LADAR.  The  difference 
being  the  APD  detectors  are  photon-counting  devices  enabling  measurements  to  be  made  at 
very  low  signal  levels  (0.4  photo-electrons  per  pixel)  as  compared  to  the  CCD  camera  (1.7 
photo-electrons  per  pixel).  The  paper  also  highlights  foliage  penetration  as  another  bene¬ 
fit  of  3D  LADAR  with  APD  detectors.  Tests  are  run  where  the  LADAR  can  see  through 
semi-transparent  material  (i.e.  camouflage  netting). 

Referring  to  [77],  a  LADAR  capability  is  presented  that  can  provide  target  informa¬ 
tion  on  sea-skimming  anti-ship  missiles.  Target  information  includes  range  and  velocity 
data.  Range  data  is  captured  by  the  time-of-flight  principle.  In  RADAR,  the  target’s  ve¬ 
locity  information  is  captured  from  the  frequency  changes  between  the  transmitted  and  re¬ 
ceived  pulses.  Typical  coherent  LADAR  architectures  require  mixing  at  light  wavelengths 
to  capture  the  differences  which  is  very  difficult.  This  paper  shows  an  interesting  work¬ 
around  combining  the  preciseness  of  laser  light  operation  and  the  mature  radio  frequency 
mixing  technology.  The  LADAR  collects  velocity  information  by  using  a  linear  frequency 
modulated  (LFM)  radio  frequency  to  amplitude  modulate  the  laser  pulse.  The  receiver  col¬ 
lects  and  coherently  mixes  in  the  RF  domain  rather  than  at  laser  light  wavelengths  thereby 
reducing  complexity. 

A  gated  3D  LADAR  is  described  where  the  detector  is  an  intensified  CCD  camera 
with  a  Nd:YAG  passively  Q-switched  32. KHz  pulsed  green  laser  at  532  nm  [6].  This  wave¬ 
length  provides  substantial  underwater  transmission.  However,  the  system  is  not  covert  or 
eye-safe  at  532  nm  like  the  SWIR  3D  FLASH  LADAR. 

In  order  to  perform  data  registration  and  extraction,  a  scanning  LADAR  is  teamed 
with  a  2D  digital  camera  [85].  This  paper  illustrates  an  example  of  using  an  active  and 
passive  system  to  increase  capabilities.  One  of  the  challenges  with  using  two  sensors  is 
fusing  the  data  sets  to  represent  the  information  in  a  consistent  coordinate  system. 

A  time-of-flight  (TOF)  “real-time”  3D  video  capability  using  a  3D  FLASH  LADAR 
is  described  in  [14].  This  paper  describes  the  architecture  required  which  is  very  similar 
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to  [81]  with  a  focal  plane  array  (FPA)  and  high-speed  ROIC  to  capture  the  range  data.  Con¬ 
sidering  GPS-denied,  GPS-degraded  regions  or  geolocation  improvement  goals,  a  method 
is  described  where  Global  Positioning  System  (GPS)  and  Inertial  Measurement  Unit  (IMU) 
data  are  fused  with  3D  FLASH  LADAR  data  by  Kalman  filtering  to  enable  autonomous  ve¬ 
hicle  control  (relative  navigation)  for  space  vehicle  docking  or  in-flight  jet  refueling. 

2. 6. 1.2  Applications.  Applications  of  the  3D  FLASH  LADAR  technol¬ 
ogy  include  target  identification,  rendezvous  operations,  foliage  penetration,  mapping,  and 
guidance  and  navigation.  Given  the  infancy  of  the  capability  and  the  interest  in  active  EO 
sensing,  this  list  will  expand  with  a  substantial  increase  in  performance  in  each  of  the  areas. 

A  comprehensive  overview  of  the  LADAR  topic  area  is  given  in  [78].  The  paper  de¬ 
scribes  utilizing  LADAR  data  to  build  synthetic  environments,  developing  LADAR  system 
models,  and  using  training  sets  for  algorithms  to  aid  in  target  recognition  and  weapon  ap¬ 
plications  (weapon  guidance,  aim  point  selection).  At  the  time,  the  authors  used  synthetic 
data  to  simulate  3D  FLASH  LADAR  data,  but  will  have  the  hardware  available  in  the  fu¬ 
ture  for  collects.  The  fusion  of  LADAR  data  with  other  sensors  yielding  impressive  results. 
Among  the  many  benefits,  one  of  the  most  important  benefits  is  more  precise  targeting 
thereby  reducing  collateral  damage. 

Using  an  innovative  scannerless  Multiple-Slit  Streak  Tube  Imaging  LiDAR  (MS- 
STIL),  [22]  reports  on  LiDAR  tests  that  demonstrate  target  imaging  through  foliage  and 
other  obscurants.  Another  test  demonstrates  capability  to  image  surf  zones  to  identify  anti¬ 
landing  mines  and  other  obstacles. 

A  variety  of  3D  scanning  LADAR  applications  are  discussed  in  [17]  relating  to  the 
use  of  APDs  in  the  receiver  design.  The  performance  of  APDs  is  reported  using  different 
materials  and  at  different  wavelengths.  Applications  include:  sensor-fused  weapons,  eye- 
safe  range-finding,  and  fire-and-forget  missiles. 

Similar  to  [14],  a  very  useful  application  for  3D  FLASH  LADAR  is  for  aerial  vehicle 
navigation  in  GPS-denied  situations  [86].  Teamed  with  IMU  data,  the  3D  FLASH  LADAR 
is  capable  of  providing  autonomous  space  vehicle  navigation  or  landing  systems  on  the 
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moon  or  other  planets.  IMU  measurements  drift  over  time  due  to  the  errors  encountered 
in  integrating  many  sensor  measurements.  GPS  is  one  mitigation  technique  to  combat  this 
drift.  In  GPS-denied  or  degraded  regions,  3D  FLASH  LADAR  data  can  replace  GPS  data 
to  recalibrate  the  IMU  measurements. 

Another  example  of  applying  3D  FLASH  LADAR  data  for  autonomous  vehicle  nav¬ 
igation  focuses  on  spaceborne  rendezvous  and  capture  [40].  The  LADAR  data  benefits  this 
application  area  by  providing  an  independent  range  to  the  docking  platform  regardless  of 
the  existence  of  other  docking  sensors.  Additionally,  the  LADAR  could  provide  an  image 
of  the  docking  platform  used  to  verify  its  integrity. 

2.6.2  3D  FLASH  LADAR  Post-Processing.  The  post-processing  of  3D  LADAR 
data  (scanning  and  FLASH  systems)  includes  range  estimation,  object  retrieval,  data  reg¬ 
istration,  edge  detection,  feature  extraction,  planar  feature  detection,  multi-sensor  assisted 
navigation  and  target  identification,  multiple  return  detection,  surface  imaging,  noise  reduc¬ 
tion,  detector  response  deconvolution,  illumination  pattern  renormalization,  jitter  removal, 
super-resolution,  and  image  enhancement.  With  the  fields  of  image  processing  and  RADAR 
being  very  mature,  the  application  of  theory  to  3D  LADARs  from  both  these  fields  is,  in 
many  cases,  novel.  While  the  processing  methods  by  themselves  are  not  new,  the  applica¬ 
tion  of  these  methods  to  the  3D  LADAR  data  set  may  have  never  been  done. 

2.6.2. 1  Range  Estimation.  In  [9],  the  waveform  parameters  (target  range, 
target  amplitude,  and  pixel  bias)  of  a  3D  FLASH  LADAR  are  estimated  via  a  maximum 
likelihood  derivation.  A  Cramer-Rao  Lower  Bound  (CRLB)  on  range  estimation  is  also 
derived.  The  unknown  target  parameters  are  estimated  by  using  maximum  likelihood  anal¬ 
ysis  on  an  idealized  waveform  model  (no  pixel  coupling).  Results  show  that  centimeter 
level  range  accuracy  is  attainable.  Closed  form  solutions  for  the  CRLB  are  provided  in 
the  follow-up  work  in  [39].  Several  different  scenarios  are  investigated  including  multiple 
returns  and  distorted  return  pulses  due  to  slanted  surfaces. 
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Referring  to  [83],  an  unusual  approach  to  range  estimation  in  a  3D  scanning  LADAR 
is  employed  called  the  Viterbi  algorithm  which  is  a  maximum  likelihood  sequence  esti¬ 
mator  (MLSE).  It  is  an  intelligent  search  algorithm  that  picks  the  most  likely  sequence  at 
each  stage  resulting  in  the  Viterbi  path.  The  Viterbi  path  results  in  an  estimated  object  from 
3D  LADAR  scans.  Without  modifying  the  algorithm,  computational  complexity  for  a  real- 
world  array  (e.g.  128  x  128)  make  this  algorithm  prohibitive.  Results  from  a  modified  VA 
algorithm  are  compared  to  a  peak  detector  and  Wiener  filter  method  showing  that  VA  out¬ 
performs  the  other  methods  in  terms  of  range  error.  The  modification  reduces  the  required 
computations. 

3D  surfaces  are  able  to  be  characterized  by  a  LADAR  system  capable  of  handling 
multiple  returns  in  a  single  received  signal  [32] .  The  LADAR  can  measure  range  and  obtain 
information  about  3D  structures  at  ranges  from  a  few  meters  to  several  kilometers.  The 
authors  employ  a  Bayesian  statistical  approach  based  on  reverse  jump  Markov  chain  Monte 
Carlo  (RJMCMC)  techniques  to  estimate  the  number,  positions  and  amplitude  of  received 
signals.  Two  types  of  receivers  are  considered  for  ranging  and  depth  measurement.  The 
types  are  Time-Correlated  Single  Photon  Counting  (TCSPC)  and  Burst  Illumination  Laser 
(BIL)  (e.g.  range  gating  or  repeated  BIL).  The  analysis  assumes  a  simplified  case  whereby 
each  pixel  is  independent  from  other  pixels.  A  Bayesian  approach  is  employed  because  it 
accounts  for  uncertainties  in  the  model  and  parameter  values  and  it  can  incorporate  prior 
knowledge  if  applicable.  A  modified  version  of  RJMCMC  incorporates  a  delayed  rejection 
step  permitting  the  Markov  chain  to  mix  better  through  different  proposal  distributions. 

Based  on  their  previous  work  [32],  the  authors  modify  an  assumptions  by  changing 
the  single  independent  pixel  model  to  one  that  includes  pixel  spatial  interdependencies  [33]. 
The  inter-pixel  dependencies  are  asserted  to  arise  naturally  in  imaging  real  world  objects. 
Again,  the  number,  positions  and  amplitudes  of  the  received  signals  are  estimated  using 
RJMCMC  incorporating  either  spatial  mode  jumping  (change  position  of  peak)  or  spatial 
birth/death  process  (creating  a  new  peak,  or  removing  a  peak). 
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Two-dimensional  range  images  are  used  to  estimate  the  target  location  and  range  [36]. 
These  estimates  are  attained  by  utilizing  a  three-dimensional  distortion  tolerant  filter  on  a 
three-dimensional  binary  representation  of  the  2D  range  image.  The  distortion  tolerant 
filter  is  derived  by  neglecting  out-of-family  correlations  and  minimizing  the  output  energy 
of  the  input  scene  due  to  additive  noise.  The  filter  is  considered  distortion  tolerant  by 
using  a  reference  target  training  data  set  to  recognize  the  targets  from  various  perspectives. 
In  [35],  the  3D  distortion  tolerant  filter  work  is  extended  to  include  the  effects  of  disjoint 
background  noise. 

In  [58],  the  authors  describe  a  3D  FLASH  LADAR  sensor  architecture  development 
with  theoretical  development  centered  around  range  processing  and  polarization  discrimi¬ 
nation  with  associated  experimental  results  attaining  range  resolutions  of  1  inch  range  res¬ 
olution  for  occluded  targets  and  0.3  inches  for  non-occluded  targets.  The  ranging  algorithm 
is  called  “bin-balancing  matched  filter”  or  BBMF™  which  uses  the  known  pulse  shape  to 
find  the  range  at  which  there  is  max  correlation  with  the  received  pulse.  A  weakness  of 
this  algorithm  is  assuming  the  transmitted  and  received  pulse  shapes  are  matched.  Sloped 
targets  and  range  clipping  makes  this  assumption  less  valid. 

The  authors  in  this  paper  use  coherent  detection  LADAR  data  with  the  expectation- 
maximization  algorithm  to  develop  a  method  to  fit  a  multi-resolution  (wavelet)  basis  to 
LADAR  range  data  in  a  maximum  likelihood  sense  [26].  The  Haar-wavelet  basis  is  used 
resulting  in  a  computationally  efficient  and  robust  algorithm.  The  wavelet  basis  is  used  for 
range  anomaly  suppression  to  decrease  range  error. 

Referring  to  [4],  a  laser  scanning  LADAR  and  several  ranging  methods  are  described. 
These  methods  include:  thresholding,  bump-hunting,  maximum  likelihood  (ML),  and  Re¬ 
versible  Jump  Markov  Chain  Monte  Carlo  Processing  (RJMCMC).  Bump-hunting  and  ML 
was  found  to  be  able  to  discern  multiple  targets  from  an  apparent  single  return.  During  low 
light  levels,  RJMCMC  was  shown  to  be  the  best  performer  in  terms  of  range  accuracy. 

2.62.2  Other  Processing  Methods.  Integration  methods  are  described 
where  3D  FLASH  LADAR  technology  is  integrated  with  inertial  measurements  to  deter- 
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mine  position  and  attitude  of  UAVs  whether  GPS  is  available  or  not  [27].  The  LADAR  data 
is  used  to  extract  planar,  line,  or  point  features  corresponding  to  walls  or  corners.  These 
features  are  combined  with  IMU  measurements  to  change  platform  attitude  or  velocity. 

In  [80],  a  3D  scanning  LADAR  is  used  to  show  the  capability  of  3D  FLASH  LADAR’s 
penetration  through  camouflage  and  foliage.  (The  authors  did  not  have  access  to  a  3D 
FLASH  LADAR  at  the  time.)  Waveform  analysis  is  performed  to  show  the  multiple  re¬ 
turn  detection  capability  important  in  FOPEN  (FOliage  PENetration).  The  Expectation- 
Maximization  algorithm  is  used  to  detect  the  number  of  peaks  in  a  given  returned  signal. 
With  the  returned  signal  assumed  to  be  a  sum  of  Gaussians,  the  mean  (target  location) 
and  standard  deviation  were  estimated.  By  using  waveform  processing,  algorithms  are  de¬ 
scribed  that  exploit  the  multiple  returns  when  the  LADAR  illuminates  vegetation  or  cam¬ 
ouflage.  By  deconvolution,  hidden  targets  under  obscuration  are  capable  of  being  detected. 
Estimation  of  target  location  and  waveform  width  is  performed  assuming  a  Gaussian  pulse 
in  a  noiseless  system,  but  no  detail  was  provided  as  to  the  estimation  method.  The  ability  to 
see  inside  a  dark  van  and  buildings  through  Venetian  blinds  is  shown.  Vegetation  removal 
to  aid  in  FOPEN  is  considered  a  research  priority  for  future  work. 

Using  the  AFRL  3D  FLASH  LADAR,  an  object  retrieval  algorithm  is  developed  for  a 
3D  FLASH  LADAR  system  illuminating  a  bar  target  using  a  microscanning  technique  [1]. 
Microscanning  is  required  in  this  system  due  undersampling  in  the  spatial  domain.  The 
microscanning  method  forces  the  eventual  data  output  to  meet  Nyquist  sampling  require¬ 
ments  by  developing  a  super  lattice  of  points  similar  to  super-resolution  techniques.  The 
object  retrieval  algorithm  was  derived  by  maximizing  the  log-likelihood  function  with  re¬ 
spect  to  a  particular  point  in  the  remote  scene  (object)  with  the  final  result  similar  to  the 
Richardson-Lucy  algorithm.  Cube  registration  (CR)  is  performed  by  computing  the  trans¬ 
lational  shifts  between  the  data  cubes  in  all  three  dimensions.  Because  the  data  cubes  are 
sampled  properly  in  the  range  dimension,  cubes  are  shifted  in  the  range  dimension  so  that 
each  cube  represents  a  common  range  to  the  target.  The  average  range  to  the  target  in  the 
data  cubes  are  calculated  and  then  compared  to  produce  a  range  so  that  each  image  frame 
within  the  data  cube  corresponds  to  the  same  distance.  In  the  spatial  domain,  transverse 
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shifts  between  cubes  are  accomplished  by  the  vector  projection  method  by  calculating  the 
global  shifts  between  corresponding  frames  in  each  data  cube  and  then  averaging  the  shifts 
across  all  frames  in  the  cube.  This  averaged  global  shift  is  assumed  the  only  shift  for  the 
cube  considering  the  fast  acquisition  time  of  the  sensor. 

In  [72],  a  scanning  laser  and  passive  electro-optical  (EO)  camera  are  used  to  create 
data  sets  enabling  sophisticated  data-processing  methods  to  use  for  building  3D  environ¬ 
ments,  data  classification,  bare  earth  extraction,  3D-reconstruction  of  buildings,  and  identi¬ 
fication  of  single  trees  and  estimation  of  their  position,  height,  canopy  size  and  species. 

Processing  methods  are  presented  that  convert  raw  3D  FLASH  LADAR  data  to  cleaned 
3D  data  cubes  enabling  information  to  be  generated,  displayed,  and  analyzed  in  real  time  [12]. 
The  processing  methods  include:  “noise  reduction,  ground  plane  identification,  detector  re¬ 
sponse  deconvolution  and  illumination  pattern  re-normalization.”  Of  most  interest  in  this 
paper  is  the  development  of  the  APD  response  deconvolution.  Ideally,  each  voxel  would 
represent  a  single  area  of  the  remote  scene.  However,  the  APD  detectors  are  not  ideal  and 
the  voxel  experience  coupling  between  each  other.  Since  the  true  APD  detector  response 
is  tough  to  measure,  the  effects  of  the  multiple-pixel  coupling  are  mitigated  by  identifying 
the  range  tails  within  the  array  and  moving  the  tail’s  energy  closer  to  the  voxel  of  interest. 

In  [10],  3D  FLASH  LADAR  data  is  used  to  collect  lacunarity  metrics  which  are  used 
to  measure  and  characterize  forest  canopy  gaps.  The  goal  is  to  establish  the  availability 
of  sub-canopy  collections  and  to  characterize  the  imaging  performance  of  different  canopy 
and  forest  types. 

Using  a  range-gated  3D  LADAR,  [79]  describes  the  ability  to  to  process  the  data  and 
characterize  different  targets  such  as  forests,  snow,  human  faces,  and  the  ability  to  penetrate 
vegetation. 

A  Bayesian  estimator  is  derived  to  perform  deconvolution  for  object  retrieval  improv¬ 
ing  3D  FLASH  LADAR  system  range  resolution  and  probability  of  detection  [7].  From  the 
deconvolution,  the  system  improves  its  ability  to  identify  surfaces  where  the  return  pulse  re- 
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fleeted  thereby  removing  the  range  estimate  ambiguity  caused  by  the  waveform  pulsewidth. 
Of  note,  no  form  of  the  object  is  specified  (i.e.  transmitted  pulse  shape). 

In  [53],  target  detection  is  performed  on  3D  LADAR  data  by  using  a  novel  3D  vol¬ 
ume  correlation  filter.  The  filter  operates  by  finding  the  parameter  value  that  maximizes 
the  volume  correlation  between  the  data  and  either  a  3D  model  or  known  3D  reference. 
Methods  of  perspective  correction  are  also  described  such  that  objects  are  represented  by 
their  true  relative  size. 

The  limits  of  theoretical  resolution  in  3D  LADAR  systems  are  derived  in  [41].  While 
previous  work  focused  on  coherent  detection  LADARs,  this  paper  extends  their  work  to 
derive  fundamental  resolution  limits  in  direct  detection  3D  scanning  and  FLASH  LADARs. 
The  “volume  of  resolution”  is  a  constant  metric  allowing  the  LADAR  designers  to  balance 
spatial  and  range  resolution  consistent  with  system  goals. 

Multiple  post-processing  methods  for  a  3D  FLASH  system  are  described  in  [35] 
including  matched  filtering,  coordinate  mapping,  jitter  removal,  and  registration.  Although, 
no  object  retrieval  methods  are  employed  to  improve  results. 

A  super-resolution  method  is  developed  for  3D  FLASH  LADAR  in  [68].  Perfor¬ 
mance  of  the  method  using  synthetic  and  real  targets  is  shown  to  be  better  than  upsampling 
and  interpolating  methods  by  using  the  Canny  edge  detection  algorithm  [11]. 

In  [15],  this  paper  develops  an  image  deconvolution  technique  using  regularized  in¬ 
version  followed  by  a  denoising  filter.  Inversion  refers  to  the  ill-posed  problem  of  removing 
the  blur  from  the  the  imaging  model.  The  inversion  process  can  produce  poor  results  in  the 
presence  of  noise  due  to  its  uniform  amplification  across  frequencies.  Regularized  inver¬ 
sion  (such  as  Wiener  filtering)  can  alleviate  such  problems.  Also,  assuming  there  are  similar 
patches  within  a  natural  image,  the  de-noising  filter  is  based  on  a  block-matching  and  3D 
(BM3D)  filtering  method.  This  work  extends  the  regularized  inversion,  regularized  wiener 
inversion,  and  BM3D  work  to  handle  colored  noise.  Of  note,  a  regularization  parameter 
that  is  determined  empirically  is  used  in  the  inversion  process. 
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2.6.3  Blind  Deconvolution.  As  part  of  the  research  effort,  blind  deconvolution 
techniques  will  eventually  be  employed  given  a  laboratory  or  field  test  with  a  3D  FLASH 
LADAR  remote  sensing  scenario.  A  review  of  the  pertinent  blind  deconvolution  literature 
is  appropriate  given  this  realistic  situation.  2D  passive  electro-optical  papers  usually  have 
one  object  and  many  different  blurring  functions  due  to  relatively  slow  image  acquisition 
times  with  corresponding  atmospheric  turbulence.  Whereas,  the  3D  FLASH  LADAR  blind 
deconvolution  scenario  has  many  different  objects  with  one  blurring  function  regardless  of 
atmospheric  turbulence  strength.  Each  data  cube  of  3D  FLASH  LADAR  is  considered  to  be 
blurred  by  one  point  spread  function  due  to  the  rapid  acquisition  time  for  the  range  images. 

Overall,  there  were  no  papers  found  that  attempted  to  restore  the  object  by  perform¬ 
ing  blind  deconvolution  on  any  type  of  3D  LADAR  system.  The  optical  astronomy  field 
dominates  the  image  blind  deconvolution  publications.  The  main  difference  between  this 
research  versus  the  typical  blind  deconvolution  is  that  this  research  endeavors  to  estimate 
the  waveform  parameters  located  within  the  object  and  a  single  point  spread  function  while 
the  typical  2D  image  blind  deconvolution  problem  estimates  the  phase  within  the  point 
spread  function  and  a  single  object.  In  other  words,  rather  than  parameterize  the  point 
spread  function,  the  object  has  been  parameterized  in  this  research. 

Generally  regarded  as  one  of  the  founding  blind  deconvolution  papers,  [59]  performs 
signal  recovery  for  multiplied  and  convolved  signals  by  using  homomorphic  filtering  uti¬ 
lizing  the  complex  cepstrum  of  the  signals.  Results  of  the  filtering  technique  applied  to 
deconvolution  problems  are  shown  in  speech  processing  and  echo  removal. 

The  other  founding  paper  concerning  blind  deconvolution  recovers  the  original  music 
from  old-time  vinyl  records  by  homomorphic  filtering  or  power  spectrum  estimation  tech¬ 
niques  [82].  The  assumed  mathematical  model  is  audible  music  resulting  from  the  original 
music  convolved  with  a  record  players  impulse  response.  They  subsequently  extend  the 
theory  to  a  simple  imaging  example  whereby  they  look  to  remove  the  effects  of  image  blur 
caused  by  camera  motion,  inaccurately  focused  lenses,  and  atmospheric  turbulence. 
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In  [42],  general  blind  deconvolution  methods  are  reviewed  and  classified  into  2  classes 
which  are  (1)  PSF  estimation  separate  from  the  true  image  estimation  and  (2)  simultaneous 
estimation  of  the  PSF  and  true  image.  The  first  class  uses  a  simple  technique  called  A  Priori 
Blur  Identification  methods.  The  second  class  incorporates  several  techniques  including 
Zero  Sheets  Separation,  Autoregressive  Moving  average  (ARMA)  Parameter  Estimation, 
Nonparametric  Deterministic  Image  Constraints  Restoration,  and  Nonparametric  Methods 
based  on  High  order  Statistics.  In  the  follow-up  paper  [43],  the  authors  discuss  other  blind 
image  deconvolution  methods  that  were  omitted  from  their  previous  article  which  were 
projection-based  blind  deconvolution  and  maximum  likelihood  restoration. 

Given  the  mathematical  model  in  this  research,  the  most  germane  article  is  from  [71]. 

This  paper  develops  a  maximum-likelihood  based  blind  deconvolution  technique  on  images 
corrupted  by  photon  noise  without  the  need  for  a  nearby  reference  point  source  which  can 
converge  to  the  solution  faster  (e.g.  less  required  frames)  than  techniques  that  do  require 
a  point  source.  The  blind  deconvolution  technique  is  called  the  Generalized  Expectation 
Maximization  (GEM)  algorithm  based  on  the  seminal  work  by  Dempster,  et  al.  [16].  The 
GEM  algorithm  is  advantageous  due  to  its  ability  to  reduce  the  maximization  complexity 
and  to  uncouple  the  object  and  blurring  function. 

In  [46],  the  blind  deconvolution  is  performed  by  error  minimization  via  conjugate 
gradient  minimization  where  the  error  is  a  composite  of  deviations  from  image  and  Fourier 
space  constraints.  Also,  blind  deconvolution  techniques  are  used  with  phase  estimation 
methods  for  object  retrieval  on  raw  speckle  images. 

Using  Kolmogorov  statistics  to  model  the  turbulent  atmosphere,  blind  deconvolution 
is  performed  on  astronomical  speckle  images  approximating  the  shot  noise  by  a  weighted 
Gaussian  noise  model  [47].  The  weighted  Gaussian  model  is  used  because  the  author  as¬ 
serts  that  many  imaging  situations  don’t  fit  the  usual  Poisson  or  Gaussian  noise  statistics. 

In  [21],  an  iterative  blind  deconvolution  algorithm  based  on  the  Richardson-Lucy  al¬ 
gorithm  is  developed  and  compared  with  a  Wiener  filter  blind  deconvolution  algorithm  [50],  [62]. 
The  authors  choose  to  develop  a  new  algorithm  based  on  the  Richardson-Lucy  algorithm 
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due  to  its  proven  robustness  in  the  presence  of  high  noise  levels.  They  also  perform  a 
“semiblind”  deconvolution  by  attempting  improve  the  algorithm  by  adding  a  priori  infor¬ 
mation  by  assuming  a  functional  form  of  the  PSF.  By  “parameterizing”  the  PSF,  the  num¬ 
ber  of  unknowns  of  the  PSF  reduces  drastically.  Conclusions  from  this  paper  are  that  the 
Richardson-Lucy  algorithm  is  more  stable  than  other  blind  deconvolution  algorithms  and 
has  a  better  noise  tolerance  than  the  Ayers-Dainty  and  Wiener  filter  algorithms.  From  [2], 
the  Ayers-Dainty  algorithm  generalizes  the  Feinup  phase  retrieval  algorithm  by  implement¬ 
ing  an  iterative  technique  based  on  Fourier  transforms  along  with  energy  conservation,  an 
image  non-negativity  constraint  and  Fourier  domain  constraints  to  estimate  the  object  and 
PSF. 

Another  attempt  to  retrieve  the  object  and  PSF  is  accomplished  by  a  maximum  a 
posteriori  (MAP)  estimator  on  a  2D  LADAR  imaging  system  [52].  Although,  in  this  case, 
it  is  the  optical  transfer  function  (OTF)  that  is  estimated  by  parameterizing  the  OTF  based 
on  Fried’s  seeing  parameter.  This  paper  also  develops  a  MAP  estimator  for  the  speckle 
parameter  in  a  negative  binomial  probability  distribution  modelling  partially  coherent  light. 

Considering  the  field  of  fluorescence  microscopy,  blind  deconvolution  is  performed 
using  an  iterative  expectation-maximization  approach  with  some  prior  knowledge  of  the 
PSF  characteristics  and  assuming  Poisson  noise  statistics  [34].  The  characteristics  include 
circular  symmetry  (general  symmetry  is  also  presented)  and  a  band-limited  nature.  The 
symmetry  argument  is  appropriate  due  to  the  symmetrical  nature  of  most  apertures.  The 
band-limit  constraint,  which  rules  out  the  trivial  solution,  also  is  appropriate  due  to  the 
low-pass  filtering  effect  of  optical  systems.  The  trivial  solution  is  the  solution  where  an  im¬ 
pulse  is  convolved  with  the  degraded  image.  Using  these  constraints,  the  algorithm  suitably 
reconstructs  the  original  images. 

In  [87],  image  recovery  is  performed  from  noisy  and  blurred  observations  by  imple¬ 
menting  an  adaptive  finite  impulse  response  filter.  Coefficients  of  this  filter  are  updated 
using  a  two-dimensional  Constant-Modulus  (CM)  cost  function  similar  to  one-dimensional 
blind  adaptive  equalization  found  in  the  communications  field. 
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Focusing  on  astronomical  applications,  this  paper  builds  on  the  iterative  blind  decon¬ 
volution  result  from  Ayers  and  Dainty  [2]  by  utilizing  methods  that  reduce  edge  effects, 
account  for  different  convergence  rates  of  the  object  and  impulse  response,  shorten  conver¬ 
gence  time,  and  perform  noise  dampening  [3].  The  methods  are  valid  when  only  constrain¬ 
ing  the  data  to  be  positive.  A  method  of  initializing  the  spatial  impulse  response  is  attained 
by  using  autocorrelations  of  the  observed  image. 

Referring  to  [18],  image  reconstruction  of  a  blurred  and  noisy  optical  system  is  per¬ 
formed  using  phase  diversity,  deconvolution  (Richardson-Lucy  based),  and  iterative  blind 
deconvolution.  All  three  methods  satisfactorily  reconstruct  the  image  with  similar  accu¬ 
racy,  but  deconvolution  is  fastest.  Their  work  handles  extended  scenes  or  scenes  in  which 
the  object  either  encompasses  the  FOV  entirely  or  is  larger  than  the  FOV.  Consequently,  the 
edge  effects  cannot  be  ignored  and  must  be  accounted  for  in  the  algorithms. 

In  [61],  blind  object  reconstruction  is  accomplished  by  reducing  the  3D  problem 
into  a  set  of  2D  problems.  Along  with  imposing  positivity  and  bandlimit  constraints,  new 
estimates  of  the  2D  image  and  PSF  are  obtained  by  Wiener  filtering  the  Fourier  transform 
of  the  image  or  PSF  respectively  with  the  current  estimate.  There  is  an  important  result 
concerning  3D  vs.  2D  sampling  requirements.  As  opposed  to  the  2D  image  scenario,  3D 
blind  deconvolution  has  a  unique  solution  even  if  the  data  is  not  Nyquist  sampled. 

2.6.4  CRB  and  Parameter  Optimization.  Compared  to  the  convolution  model 
contained  in  the  present  paper,  previous  work  on  bounding  range  performance  in  the  LADAR 
topic  area  focused  on  single  pixel  (i.e.  single  target  in  a  pixel)  analysis.  In  [9],  a  CRB  on 
range  estimation  is  derived  for  a  single  pixel  of  a  3D  FLASH  LADAR.  In  support  of  the 
bound,  the  unknown  waveform  parameters  (target  range,  target  amplitude,  and  pixel  bias) 
are  estimated  via  a  maximum  likelihood  estimation  algorithm.  Theoretical  and  simulation 
results  show  that  centimeter  level  range  accuracy  is  attainable.  Closed  form  solutions  for 
the  CRB  are  provided  in  the  follow-up  work  in  [39]. 

Another  paper  developed  a  signal-to-noise  (SNR)  based  method  to  determine  range 
and  spatial  resolution  limits  of  scanning  and  direct  detection  LADAR  [41].  While  account- 
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ing  for  the  proper  LADAR  noise  sources  and  operating  parameters,  the  SNR-based  method 
does  not  consider  the  performance  of  the  algorithms  required  to  estimate  the  resolution  in 
the  presence  of  noise. 

Other  literature  has  utilized  the  Gaussian  function  to  describe  the  object.  In  [28], 
the  object  profile  is  defined  by  a  Gaussian  in  one  dimension  corrupted  by  additive  Gaussian 
noise.  The  CRB  on  a  one  target  profile  estimation  is  performed.  In  another  paper,  the  object 
is  a  two-dimensional  (2D)  Gaussian  describing  the  incident  intensity  on  a  charge-coupled 
device  (CCD)  array  [89].  This  2D  Gaussian  is  used  to  develop  a  two-dimensional  CRB 
for  any  unbiased  position  estimator  as  well  as  a  maximum-likelihood  (ML)  position  optical 
estimator  (position  only,  no  range  information  or  estimate). 

The  use  of  the  CRB  in  parameter  optimization  or  performance  characterization  has 
been  done  previously  in  fields  such  as  heterodyne  Light  Detection  and  Ranging  (LiDAR), 
RADAR,  and  positron  emisson  tomography  (PET)  [20],  [44],  [48],  [49],  and  [69].  In  all  the 
papers,  the  method  was  to  pick  the  optimum  condition  based  on  CRB  minimization  either 
through  physically-based  analytic  expressions  or  bound  comparisons  over  different  parame¬ 
ter  choices.  In  [69],  comparisons  are  made  using  the  CRB  concerning  Doppler  estimation  in 
heterodyne  and  direct  detection  LiDAR  given  several  different  operating  parameters.  Also, 
methods  are  discussed  that  enable  heterodyne  Doppler  estimation  performance  to  approach 
that  of  the  CRB.  Concerning  synthetic  aperture  RADAR  (SAR)  design  in  [49],  the  CRB 
developed  in  this  paper  showed  that  performance  is  enhanced  by  implementing  a  multi¬ 
dimensional  aperture  over  a  one-dimensional  aperture.  In  [48],  the  CRB  is  used  to  validate 
the  use  of  range  compression  in  multi-input  multi-output  (MIMO)  RADAR.  Also,  wave¬ 
form  optimization  in  MIMO  RADAR  is  accomplished  via  several  minimization  techniques 
on  the  CRB  matrix  to  include  minimizing  the  trace,  determinant  and  largest  eigenvalue. 
Another  paper  uses  the  CRB  to  select  an  optimal  RADAR  beamspace  transformation  oper¬ 
ator  [20].  The  optimality  condition  metric  is  physically-based  using  the  analytical  form  for 
the  beamspace  transformation  that  minimizes  the  CRB  function  itself.  Finally,  the  design 
parameters  of  avalanche  photo  diodes  (APD)  used  in  small  animal  PET  are  optimized  by 
selecting  those  parameters  from  the  search  space  that  the  minimize  the  CRB  [44]. 
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III.  Laboratory  Data  Collection 

In  order  to  verify  theory  and  simulation  range  estimation  results,  laboratory  measure¬ 
ments  were  collected  using  an  Advanced  Scientific  Concepts  (ASC)  Inc.  three  dimen¬ 
sional  FLASH  LAser  Detection  And  Ranging  (3D  FLASH  LADAR)  that  illuminated  a 
target  corresponding  to  one  used  in  simulation.  The  details  behind  the  collection  are  the 
topic  of  this  chapter. 

Using  the  three  bar  target  template,  a  laboratory  experiment  was  conducted  using  3D 
FLASH  LADAR  hardware  consistent  with  parameters  in  Table  3.1.  Experimental  results 
presented  in  a  later  chapter  show  range  estimation  improvement  after  applying  the  object 
recovery  techniques.  However,  modifications  to  the  camera  and  raw  data  were  necessary  to 
enable  a  proper  experiment  and  ensure  that  the  data  matches  the  model  from  Section  2.5. 
The  system  point-spread-function  (PSF)  is  also  determined  experimentally  using  a  step 
target  which  is  done  such  that  the  PSF  can  be  used  in  the  object  decoloration  algorithm 
(Wiener  filter)  detailed  in  Section  5.1.1.  Finally,  the  ability  to  use  object  recovery  algo¬ 
rithms  is  contingent  on  using  the  incoherent  light  model  described  in  Section  2.1.  Thus, 
the  speckle  parameter  of  the  partially  coherent  light  distribution  is  estimated  and  compared 
against  the  incoherent  model.  While  some  speckle  noise  is  evident  in  the  data,  the  estima¬ 
tion  results  indicate  that  the  incoherent  model  is  a  valid  approximation. 

The  chapter  is  organized  as  follows:  Section  3.1  provides  details  on  the  3D  FLASH 
LADAR  hardware,  Section  3.2  discusses  the  laboratory  collection  set-up  used  for  experi¬ 
mental  data  processing  in  Section  5.3,  Section  3.3  identifies  the  default  hardware  configura¬ 
tion  as  spatial  aliased  and  describes  the  correction,  Section  3.4  provides  the  steps  required  to 
pre-process  the  experimental  data  including  gain  variation  equalization  and  photon  scaling, 
Section  3.5  specifies  how  the  system  PSF  was  attained  from  a  step  target,  and  Section  3.6 
derives  a  speckle  parameter  estimator  and  performs  the  estimation  on  the  experimental  data. 

3.1  3D  FLASH  LADAR  Hardware  Description 

A  3D  FLASH  LADAR  is  an  active,  pulsed  system  that  is  both  an  imaging  and  rang¬ 
ing  sensor.  It  produces  a  time  sequence  of  two-dimensional  (2D)  images  due  to  a  fast  range 
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Table  3.1:  3D  FLASH  LADAR  parameters 


Parameter 

Value 

Detector  Array 

128  x  128 

Aperture  Diameter  (D) 

2  mm 

Mean  Wavelength 

1.55  //m 

Focal  Length 

0.30  m 

Target  Range 

5.21  m 

Transmit  Energy 

10  mJ 

Pulse  Standard  Deviation  (aw) 

3  ns 

Beam  Divergence 

0.009  radians 

Detector  Spacing 

100  jum 

Detector  Array  Fill  Factor 

100% 

Detector  Bandwidth 

0.5  /im 

Target  Reflectivity 

10% 

Solar  Irradiance 

10  Watts/ni2///m 

D/ra  Seeing  Condition 

1.43 

Frame  Rate 

30  Hz 

Time  Samples 

20 

Sample  Period 

1.876  ns 

gate  resulting  in  a  3D  data  cube  of  spatial  and  range  scene  data  with  excellent  range  reso¬ 
lution  [19],  [81].  FLASH  technology  principally  differs  from  scanning  LADAR  by  being 
able  to  form  a  3D  representation  of  a  remote  scene  in  one  laser  pulse  rather  than  rastering 
a  3D  scene  together  using  many  pulses.  This  capability  results  in  faster  scene  collection 
times  with  lighter  weight,  lower  power,  and  reduced  mechanical  complexity  as  compared 
to  the  scanning  systems.  Operating  in  the  eye-safe  short-wave  infrared  region  (SWIR)  of 
the  electromagnetic  spectrum  at  1570  nm,  a  representative  system  shown  in  Figure  3.1  is 
built  by  ASC,  Inc.  and  has  receiver  electronics  consisting  of  a  128  x  128  detector  array  and 
associated  circuity  capable  of  producing  twenty  (20)  2D  range  slice  images  [66].  Detector 
pixel  separation  is  100  micrometers  with  nearly  100%  fill  factor  due  to  a  focusing  micro¬ 
lens  array  in  front  of  the  detector  pixel  array.  An  extremely  fast  range  sampling  interval 
of  1 .876  ns  makes  it  nearly  impervious  to  platform  motion  distortion  for  a  single  cube  col¬ 
lection.  Depending  on  the  LADAR  operating  mode,  each  pixel  could  either  have  a  distinct 
starting  range  dependent  on  received  photon  levels  (“hit  mode”)  or  have  the  same  starting 
range  (“sular  mode”).  Capable  of  “real-time”  3D  movies,  it  produces  a  cube  of  spatial  and 
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Figure  3.1:  A  picture  of  the  Applied  Scientific  Concepts  (ASC)  Inc.  3D  FLASH  LADAR 
system  including  the  laser,  receiver  optics  and  electronics,  and  laptop.  ASC 
provides  a  laptop  to  operate  the  LADAR  and  view  and  process  the  received 
signals. 


range  scene  data  where  each  2D  range  slice  image  contains  the  detected  counts  proportional 
to  the  incident  photo-electrons  upon  each  pixel  in  the  detector  array.  Four  dimensions  of 
data  are  available  to  include  the  two  spatial  coordinates,  range,  and  intensity. 

As  previously  described,  a  3D  FLASH  LADAR  operates  in  one  of  two  modes.  The 
first  mode  is  called  “hit  mode”  where  each  pixel  element  (pixel)  is  independently  triggered 
when  its  intensity  reaches  a  preset  threshold.  This  mode  is  advantageous  when  searching  for 
a  target  where  the  range  is  not  already  known.  However,  truncated  waveforms  can  occur 
leading  to  range  estimation  errors.  The  second  mode  is  called  “sular  mode”  where  the 
pixels  are  triggered  to  start  recording  data  together  based  on  a  preset  range.  Benefits  of  this 
mode  include  being  able  to  successively  capture  fine  details  of  the  target  and  background. 
Drawbacks  are  that  the  target  range  must  be  known  a  priori  and  waveforms  are  truncated 
for  targets  near  the  end  of  the  collect.  An  potential  CONcept  of  Operations  (CONOP)  is  for 
“hit  mode”  to  operate  like  a  search  RADAR  and,  once  the  target  is  acquired,  “sular  mode” 
would  track  and  identify  the  target. 
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The  breakthrough  technology  in  the  ASC  3D  FLASH  LADAR  is  the  Laser  RADAR 
Processor  (LRP)  which  allowed  for  the  fast  range  sampling  and  independent  pixel  con¬ 
trol  [63].  Due  to  advances  in  semiconductor  technology,  the  LRP  was  originally  a  32  x 
32  detector  array  with  400  pm  pixel  separation  which  improved  to  a  128  x  128  array  with 
100  pm  pixel  separation  using  Indium  Gallium  Arsenide  (InGaAs)  avalanche  photo-diodes 
(APD)  as  the  detector  material.  APD  detectors  generate  many  electrons  for  a  single  inci¬ 
dent  photon  and  are  useful  in  low-light  situations.  The  fast  range  sampling  is  achieved  by 
analog  and  digital  circuitry  independently  located  behind  each  of  the  pixels. 

3.2  Data  Collection  Details 

Located  at  Wright-Patterson  AFB,  OH,  the  Air  Force  Research  Laboratory  (AFRL) 
Sensors  Directorate  contains  facilities  acceptable  for  operation  of  the  3D  FLASH  LADAR. 
Ideally,  the  intent  would  be  to  operate  the  LADAR  from  the  top  floor  of  the  building  across 
a  considerable  distance  (kilometers)  right  after  dusk  to  experience  atmospheric  turbulence. 
However,  due  to  the  constraints  of  the  aperture  size  (discussed  in  Section  3.3),  the  target 
range  is  shortened  to  meters  to  allow  for  a  sufficient  signal-to-noise  ratio  (SNR).  The  range 
to  the  first  surface  is  5.21  m  and  is  set  up  to  be  1.7  m  into  the  range  collections.  Range 
to  the  second  surface  is  1 .22  m  from  the  first  surface  to  give  roughly  four  range  samples 
between  surfaces. 

Receiver  optics  required  some  modifications  from  the  default  configuration  [66],  [9]. 
The  optics  are  focused  on  the  first  surface  which  means  that  the  successive  range  collects 
are  slightly  de-focused.  The  resulting  data  shows  little  effect  from  the  lack  of  focus.  Con¬ 
sidering  the  short  range  distance,  a  one  degree  diffuser  is  put  on  the  laser  transmission 
optics  to  enable  the  entire  target  to  be  illuminated  by  the  beam  without  lowering  the  SNR 
prohibitively.  The  focal  length  is  set  at  300  mm.  Due  to  sampling  issues  covered  in  a  later 
section,  the  aperture  diameter  is  changed  to  2  mm  by  using  brown  cardboard  with  a  circu¬ 
lar  hole  cut  in  the  center  placed  in  front  of  the  supplied  aperture  (10  cm).  Using  similar 
triangles,  each  detector  pixel  field  of  view  (FOY)  xp  corresponds  to  1.7  mm  at  the  target 
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Figure  3.2:  True  ranges  of  the  three-bar  target  with  first  surface  at  5.21  m  and  second 
surface  at  6.43  m  with  1.22  m  of  separation  in  between  surfaces. 

location  determined  by  the  following  calculation: 

Xp  Xd 

Xt  fl 

Xp  100  n m 

5.21  m  300  mm 

Xp  =  1.7  mm  (3.1) 

where  xt  is  the  target  range,  x,i  is  the  pixel  separation,  and  fi  is  the  focal  length. 

Referring  to  target  template  depicted  in  Figure  3.2,  the  first  surface  targets  are  con¬ 
structed  from  white,  flat  cardboard  with  the  bars  cut  out  of  one  board  (first  surface)  and  the 
other  board  is  left  untouched  (second  surface).  There  are  two  slimmer  rectangle  targets  and 
one  larger  rectangle  target.  The  slimmer  targets  are  0.5  cm  width  by  5  cm  length  and  the 
larger  target  is  1  cm  width  by  5  cm  length.  All  three  targets  are  individually  separated  by 
1.5  cm  (edge  to  edge). 

3.3  Spatial  Aliasing 

Due  to  limits  in  current  detector  technology  requiring  a  large  footprint  for  the  elec¬ 
tronics  behind  each  pixel,  the  receiver  optics  are  spatially  under-sampled  which  needs  to  be 
mitigated  in  order  for  the  received  data  to  be  unaliased.  The  aliasing  would  cause  uncer¬ 
tainty  in  the  received  data  and  violate  the  data  model.  The  default  configuration  is  aliased 
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because  of  the  Nyquist  sampling  theory  in  which  the  sampling  rate  must  be  at  least  twice 
the  highest  frequency  content  in  the  signal.  The  optics  are  a  natural  low-pass  filter  with 
the  highest  frequency  called  the  cut-off  frequency.  For  incoherent  imaging,  the  cut-off 
frequency  is  [25] 

fo  =  TT  (3-2) 

AZi 

where  D  is  the  aperture  (exit  pupil)  diameter,  A  is  the  light  wavelength,  and  zt  is  the  image 
distance.  Therefore,  the  focal  plane  must  sample  at  twice  this  spatial  frequency  or  2D/\Zi. 
The  typical  apertures  for  this  camera  are  in  the  centimeters.  For  example,  an  aperture  of  10 
cm  would  equate  to  a  spatial  frequency  sampling  requirement  at  4.3xl05  cycles  per  meter. 
At  100  /im  spacing,  the  detector  array  does  not  meet  this  requirement.  If  the  aperture 
is  reduced  to  2  mm,  then  the  spatial  frequency  sampling  requirement  is  now  at  8.6xl03 
cycles  per  meter  which  the  detector  array  can  meet.  However,  the  aperture  reduction  comes 
at  the  expense  of  reduced  light  gathering  and  shortened  range  in  which  the  LADAR  can 
be  operated.  Thus,  the  target  range  is  placed  at  5.21  meters  (near  the  minimum  ranging 
distance  of  the  sensor)  to  obtain  high  enough  signal  to  noise  ratio  (SNR)  in  the  collected 
data. 

3.4  Data  Pre-processing 

The  data  observations  from  the  3D  FLASH  LADAR  hardware  need  pre-processing 
steps  to  be  suitable  for  insertion  into  the  Wiener  filter  and  GEM  algorithms.  In  simu¬ 
lation,  the  noisy  and  blurry  data  are  well-controlled  and  therefore,  well-behaved.  While 
the  experimental  3D  FLASH  LADAR  data  exhibits  expected  pixel  waveform  shapes  (i.e. 
Gaussian-like)  and  spatial  blur,  the  data  is  ill -behaved  to  a  degree  due  to  inherent  features 
of  the  hardware  performance. 

Referring  to  [38]  and  [73],  the  experimental  hardware  experiences  a  gain  phenomenon 
whereby  a  pixel’s  gain  drops  when  laser  energy  is  incident  upon  a  large  area  of  another  part 
of  the  detector  array.  With  the  three  bar  target,  the  laser  energy  is  incident  on  the  front 
surface  first  which  causes  second  surface  pixels  to  experience  a  gain  drop.  Figure  3.3(b) 
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Figure  3.3:  (a)  Gain  profile  correction  resulting  from  executing  Equation  (3.7).  By  look¬ 

ing  at  background  pixels,  the  hardware  gain  dip  is  clearly  evident  at  the  first 
surface  (near  range  sample  five)  and  the  second  surface  (near  range  sam¬ 
ple  nine).  The  first  surface  gain  drop  is  larger  than  the  second  surface  gain 
drop  due  to  the  larger  number  of  pixels  illuminated  (i.e.  larger  surface  area). 
Amount  of  gain  drop  is  proportional  to  received  intensity  level  and  quantity 
of  pixels  illuminated. 

(b)  Investigating  Pixel(  19,32)  from  experimental  three  bar  target,  the  pixel 
waveform  benefits  from  the  gain  variation  correction  by  removing  the  gain 
drop  near  range  sample  four.  After  correction,  the  pixel  waveform  looks 
more  like  the  intended  pulse  model,  but  with  unwanted  noise  artifacts. 
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shows  the  gain  drop  for  a  second  surface  pixel.  The  method  for  correcting  the  gain  is  to 
calculate  an  average  gain  profile  by  looking  at  background  pixels  (i.e.  returned  laser  energy 
not  incident  on  these  pixels). 

Assuming  the  system  noise  follows  the  Poisson  distribution  and  the  gain  is  constant 
between  pixels,  the  data  model  for  an  arbitrary  pixel  is 


d(t)  =  G(t)  [Is(t)  +  IB(t)} 


(3.3) 


where  G  (£)  is  the  unitless,  time- varying  gain,  Is  (t)  is  the  laser  signal  in  units  of  photons, 


and  IB  (£)  is  the  background  signal.  A  new  variable  d  (t)  is  determined  by 


(3.4) 


where  iB  (£)  =  G  (t)  I B  (t)  is  a  known  average  background  signal  with  gain  and  IB  (t)  is 
the  mean  background  signal  without  gain.  The  variable  i B  (t)  is  separately  calculated  in 
the  laboratory  by  averaging  the  detected  background  signal  for  selected  voxels  across  many 


data  cubes.  Looking  at  the  background  pixels  only,  d  (£)  is 


G(t)IB(t)  _  G(t)IB(t)  _  IB(t ) 


(3.5) 


iB  ( t )  G  (t)  IB  (t)  IB  (■ t ) 


Taking  the  statistical  variance  results  in 


(3.6) 


Applying  this  result  and  using  a  sample  variance  of  d  in  place  of  the  statistical  variance 
(s2  — y  var  ( d(t ) ) ),  the  gain  is  determined  by 


IB  (tj 


(3.7) 
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and  can  be  seen  in  Figure  3.3(a).  This  gain  profile  is  used  on  each  of  the  pixel’s  waveforms 
to  correct  for  the  hardware  deficiencies  and  to  more  closely  match  the  model.  For  example, 
Figure  3.3(b)  shows  the  benefits  of  the  gain  correction  for  one  second  surface  pixel.  Also 
observed  in  the  previous  work,  a  side  benefit  of  gain  correction  in  both  first  and  surface 
pixels  is  the  waveform  becomes  more  symmetrical.  The  emitted  laser  pulse  shape  is  a 
hybrid  of  a  Gaussian  or  negative  parabolic  shape  with  some  asymmetry.  Gain  correction 
takes  out  some  of  the  asymmetry. 

The  3D  FLASH  LADAR  is  also  not  a  photon-counting  device  where  one  digital  count 
equals  one  photon.  The  receiver  optics  use  Avalanche  Photo  Diodes  (APD)  where  one 
photon  equals  many  detected  counts.  Consequently,  intensity  scaling  must  be  performed 
to  condition  the  data  to  be  consistent  with  the  Poisson  distribution.  The  conditioning  is 
performed  by  using  the  statistics  of  the  light  and  the  detected  mean  and  variance  of  the 
data.  The  detected  mean  of  the  data  is  qK  where  q  is  a  scaling  factor  with  units  of  detected 
counts  per  photon  and  K  is  the  true  mean  in  units  of  photons.  Since  incoherent  imaging  is 
assumed,  the  detected  variance  becomes 

q2a2  =  q2  K  (3.8) 

noting  that  the  mean  and  variance  of  the  Poisson  distribution  are  the  same.  The  data  is 
scaled  by  solving  for  q  and  then  converting  the  detected  counts  to  photons  by 

dph  =  —  (3.9) 

q 

where  dph  is  the  data  in  units  of  photons  and  ddc  is  the  data  in  units  of  detected  counts. 

3.5  Experimental  PSF 

The  Wiener  filter  is  used  to  provide  a  comparison  to  the  GEM  algorithm  [55].  In  order 
to  implement  the  Wiener  filter,  the  PSF  must  be  known.  Since  the  derivative  of  a  system  step 
response  is  the  system  impulse  response,  the  PSF  is  determined  by  taking  the  derivative  of  a 
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experimental  step  target.  Figure  3.4(a)  shows  a  range  image  of  the  step  target  collected  with 
the  same  hardware  as  the  bar  target  data.  Although,  the  entire  range  image  does  not  meet 
the  requirements  of  being  a  step  target  due  to  the  non-uniform  intensity  on  the  left-hand- 
side  (LHS).  Therefore,  a  symmetric  impulse  response  was  assumed  and  the  right-hand-side 
(RHS)  of  the  impulse  response  was  copied  and  flipped  over  to  use  as  the  LHS.  Figure  3.4(b) 
exhibits  the  resulting  profile  with  an  outer  product  operation  producing  the  two-dimensional 
PSF.  Phase  retrieval  is  then  performed  via  the  Gerchberg-Saxton  algorithm  to  arrive  at  the 
PSF  used  by  the  Wiener  filter  [23].  This  requirement  to  know  the  PSF  is  a  shortcoming 
of  the  Wiener  filter  algorithm.  Figures  3.4(c)-(d)  show  the  optical  transfer  function  (OTF) 
where  the  optics  exhibit  a  nearly  diffraction-limited  performance. 

3. 6  Speckle  Parameter  Estimation  -  Incoherent  Imaging 

Both  the  negative  binomial  and  Poisson  distributions  can  be  used  to  capture  the  non¬ 
negative,  discrete  nature  of  the  laser  light.  The  negative  binomial  distribution  would  be 
the  most  optimal  in  describing  the  illuminating  partially  coherent  laser  light,  but  blind  de- 
convolution  methods  are  cumbersome  [24].  Whereas,  blind  deconvolution  methods  with 
the  Poisson  distribution  (incoherent  imaging)  are  more  tractable  and,  thus,  utilized  in  this 
research.  Even  if  the  speckle  is  severe,  the  benefit  of  modeling  the  speckle  does  not  out¬ 
weigh  the  cost  of  implementing  a  partially  coherent  blind  deconvolution  model  for  the 
3D  FLASH  LADAR  system.  Previous  research  using  the  incoherent  data  model  for  a  3D 
FLASH  LADAR  has  also  experienced  success  [9],  [39]. 

To  gain  more  insight  into  this  assumption,  a  simple  approach  is  to  estimate  the 
amount  of  coherence  contained  within  the  3D  LLASH  LADAR  data  by  estimating  the 
speckle  parameter  of  the  negative  binomial  distribution  directly  from  the  data  [24].  Captur¬ 
ing  both  temporal  and  spatial  coherence,  if  the  speckle  parameter  estimate  is  high  enough, 
the  negative  binomial  distribution  wifi  look  Poisson-like  allowing  the  data  observations 
to  be  modeled  as  arising  from  an  intensity  convolution  (incoherent  imaging).  Including 
speckle  and  photon  noise  effects,  the  negative  binomial  probability  mass  function  (PML) 
describes  the  photon  distribution  of  a  partially  coherent  imaging  system  for  a  single  pixel 
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Figure  3.4:  (a)  One  range  image  of  the  step  target  data  cube.  Although  the  board  edge 

is  clearly  visible,  the  variable  intensity  across  it  causes  an  issue  with  the 
impulse  response  calculation.  The  step  response  definition  requires  a  con¬ 
stant  amplitude  at  all  spatial  positions.  The  target  board  portion  of  the  step 
response  does  not  meet  this  requirement,  but  the  non-target  area  (right-hand- 
side)  does  exhibit  a  constant  amplitude.  The  portion  of  the  step  response 
function  where  it  turns  off  is  this  non-target  area.  Performing  the  step  re¬ 
sponse  derivative  only  on  this  non-target  area  solves  the  problem  of  variable 
target  board  amplitude. 

(b)  ID  cut-out  of  the  resulting  PSF.  Assuming  circular  symmetry,  an  outer 
product  operation  is  used  to  find  the  corresponding  2D  PSF. 

(c)  Optical  transfer  function  (OTF).  The  OTF  is  found  by  taking  the  Fourier 
Transform  of  the  experimental  PSF  [25]. 

(d)  ID  cut-out  (zero  spatial  frequency)  of  the  OTF.  The  profile  shows  nearly 
diffraction-limited  optics  with  a  cut-off  frequency  at  4050  cycles  per  meter. 
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or  [24] 


P(K) 


T(K  +  M) 

r(A^  +  i)r(2U) 


K 
+  M 


-M 


(3.10) 


where  Ai  is  the  speckle  parameter  and  K  is  the  pixel’s  average  photon  count.  Changing  the 
distribution  for  a  3D  FLASH  LADAR,  the  illuminating  laser  light  statistics  for  a  particular 
volume  element  (voxel)  (x,y,k  remain  constant)  across  many  data  cubes  is 


P  (Djk  (x,  y)  =  djk  (x,  y)  Vj  e  (1,  2, ...,  J))  = 


nr  (djk  (x,  y)  +  M) 
J=1T  (djk  (x,y)  +  l)T(M) 


„  M 
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(3.11) 


where  j  represents  the  data  cubes,  k  is  the  range  image  (i.e.  time  variable)  within  a  data 
cube,  (x,  y)  are  the  coordinates  in  the  image  plane,  and  d)k  (x.  y)  is  the  data  observation. 
The  voxels  are  assumed  statistically  independent  from  each  other  because  of  the  discrete 
nature  of  photons  and  the  detected  photons  do  not  affect  future  detected  photons.  The 
maximum  likelihood  solution  for  the  average  voxel  intensity  is  determined  by 


K 


-j^djk  (x,y). 


3= 1 


(3.12) 


Taking  the  natural  log  of  Equation  (3.11)  yields 


In  [P  ( Djk  (x,  y)  =  djk  (x,  y)  Vj  e  (1,  2, ...,  J))]  = 
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(3.13) 


where  graphical  methods  are  employed  to  find  the  speckle  parameter  that  maximizes  this 
log-likelihood.  Using  the  same  experimental  data  as  in  the  range  estimation  efforts,  a  col¬ 
lection  of  voxels  with  the  strongest  laser  light  is  chosen  to  estimate  the  speckle  parameter. 
Figure  3.5  shows  the  similarities  between  the  negative  binomial  and  Poisson  distribution 
using  an  average  of  the  estimated  speckle  parameter. 
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Figure  3.5:  This  plot  shows  the  negative  binomial  (NB)  using  an  estimated  average 

speckle  parameter  (M  =  414)  versus  the  Poisson  distribution  with  the  same 
mean  ( K  =  3447).  While  not  identical,  the  negative  binomial  distribution 
compares  well  enough  to  the  Poisson  distribution  to  assume  incoherent  imag¬ 
ing. 

Even  without  considering  speckle  parameter  estimation  results,  the  argument  can  be 
made  for  incoherent  imaging  due  to  the  Poisson  distribution’s  ability  to  model  the  non¬ 
negativity  and  discrete  nature  of  light  [9],  [39].  This  argument  is  solidified  by  the  speckle 
parameter  estimation  results  indicating  that  the  speckle  noise  appears  low  enough  for  the 
incoherent  imaging  model  to  be  used  with  confidence. 
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IV.  Range  Estimation 


Range  estimation  in  three  dimensional  FLASH  LAser  Detection  And  Ranging  (3D 
FLASH  LADAR)  has  been  limited  thus  far  to  statistical  methods  that  operated  on 
data  models  that  did  not  incorporate  the  blurring  effect  of  the  spatial  impulse  response.  In 
other  words,  there  was  a  one-to-one  mapping  between  the  object  plane  and  image  plane 
points.  Considering  a  the  3D  FLASH  LADAR  system  as  a  linear,  space-invariant  process, 
the  relationship  between  the  object  and  image  plane  is  fully  described  by  a  convolution 
between  the  object  plane  intensity  and  the  intensity  point  spread  function  (spatial  impulse 
response).  Consequently,  the  simple  models  ignore  the  pixel-to-pixel  coupling  that  could 
significantly  degrade  range  estimation  results.  Referring  to  Chapter  II,  whether  using  a 
simplified  model  or  a  higher  fidelity  model,  the  method  of  ranging  is  exactly  the  same  in 
that  each  pixel  in  the  detector  array  is  ranged  independently. 

This  chapter  details  several  pixel-based  ranging  algorithms  include:  Section  4.1  - 
peak  detection,  Section  4.2  -  maximum  likelihood  [55],  Section  4.3  -  normalized  cross¬ 
correlation  [56],  Section  4.4  -  two  point  target  estimator,  and  Section  4.5  -  two  surface 
estimator.  The  two-point  target  estimator  is  a  novel  contribution  that  is  able  to  spatially  and 
temporally  estimate  two  point  targets  in  a  scene. 

4.1  Peak  Detection 

A  very  simple  ranging  algorithm  is  peak  detection.  This  algorithm  selects  the  range 
sample  D  (x,  y )  based  on  where  the  peak  sample  count  occurs  or 

D  (x,y)  =  argmax4  (x,y) .  (4.1) 

k 

where  <4  (x,  y)  is  the  received  waveform,  k  is  the  range  sample  variable,  and  (x,  y)  are 
the  pixel  dimensions.  Theoretically,  if  the  received  waveform  was  sampled  continuously, 
one  could  perform  peak  detection  and  not  encounter  any  interpolation  or  quantization  error. 
However,  real  systems  have  a  sampling  period  which  creates  some  ambiguity  when  peak 
detection  is  used.  Therefore,  more  capable  methods  are  sought  that  enable  estimation  to 
be  sub-sample.  Some  of  the  errors  though  could  be  mitigated  by  interpolation.  The  effects 
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of  spatial  coupling  and  shot  noise  would  further  degrade  the  waveform  of  a  sub-sample 
target  in  addition  to  the  deformation  already  encountered  by  its  sub-sample  range  position. 
These  effects  would  make  obtaining  accurate  estimates  from  standard  peak  detection  very 
difficult. 


4.2  Maximum  Likelihood 

Based  on  [55],  this  section  reviews  the  development  of  a  maximum  likelihood  method 
to  estimate  range  to  the  target  at  a  single  pixel  given  a  transmitted  Gaussian  pulse  with  ad¬ 
ditive  Gaussian  noise.  Maximum  likelihood  is  chosen  because  of  its  relation  to  the  Gener¬ 
alized  Expectation  Maximization  (GEM)  algorithm  used  in  a  subsequent  chapter  where  an 
iterative  technique  possibly  leads  to  the  maximum  likelihood  solution.  From  [84],  the  max¬ 
imum  likelihood  estimator  is  the  parameter  estimate  where  the  maximum  of  the  a  posteriori 
density  occurs.  Using  Gaussian  statistics  to  describe  the  incoming  noise,  this  a  posteriori 
density  for  an  arbitrary  pixel  (x,  y )  and  range  sample  k  is 


P[Dk(x,y)  =  dk(x,y)\ 


1  -(dkix’y'i-lk(x'y))2 

. —  e  2^2 
a/27 rex 


(4.2) 


where  a  is  the  Gaussian  noise  standard  deviation.  The  remaining  derivation  assumes  de¬ 
pendence  on  (x,  y)  and  drops  the  notation.  Since  there  are  N  time  samples  and  assuming 
the  time  samples  are  statistically  independent  of  each  other,  the  total  distribution  across  all 
time  samples  is  the  product  of  the  individual  distributions: 


P[Dk 


K 

dk-\/k  e  [1,  ...,K]]  = 

k= 1 


~( dk~ik'> 2 

e  2^ 


(4.3) 


Given  that  maximizing  the  natural  log  of  a  function  is  the  same  as  maximizing  the  function 
itself,  taking  the  natural  log  of  Equation  (4.3),  L  =  In (P(d(tk)),  results  in  the  advantageous 
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form 
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Because  the  range  and  amplitude  are  both  unknown  parameters,  the  estimation  process 
must  estimate  the  amplitude  first  and  is  found  by  [84] 


dmi(R)  =  argmax  L. 

A 


(4.5) 


Taking  the  derivative  with  respect  to  A  in  Equation  (4.4)  and  setting  it  equal  to  zero  results 
in 


k=i 


Apk(R)  -  B ) 
2a2 


Pk(R)  =  0 


(4.6) 


where  the  term  that  doesn’t  depend  on  A  has  been  dropped.  Grouping  terms  and  canceling 
a2  gives 

K 


([4  -  B]pk(R)  -  Apl(R ))  =  0  (4.7) 

k=\ 


Solving  for  the  amplitude  of  the  received  waveform,  A,  results  in 


E  [dk-B]pk(R)) 

dmi(R)  =  - •  (4.8) 

EpI(R) 

k= 1 


One  important  observation  of  the  Equation  (4.8)  is  its  dependence  on  range.  For  each 
pixel,  the  amplitude  estimation  process  consists  of  selecting  a  candidate  range  R  in  pk  and 
stepping  through  each  time  sample  to  determine  the  maximum  likelihood  solution  for  A. 
Using  this  amplitude  estimate,  the  only  other  unknown  for  a  given  pixel  is  the  target  range. 
Finding  a  similar  closed  form  solution  for  a  range  estimate  is  troublesome  due  to  the  range 
term  residing  in  the  exponential.  Hence,  finding  the  maximum  of  the  distribution  with 
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respect  to  range  by 


fmi  =  arg  max  L  (4.9) 

R 

is  mathematically  equivalent  to  using  the  amplitude  estimate  to  calculate  the  values  for  L 
in  Equation  (4.4)  for  each  candidate  range  and  selecting  the  range  that  corresponds  to  the 
largest  L  value.  This  range  serves  as  the  estimated  range  for  that  pixel.  The  algorithm  for 
estimating  the  range  in  each  pixel  is  thus: 

1)  Select  pixel  location  (x,  y) 

2)  Select  candidate  range 

3)  Estimate  waveform  amplitude,  A 

4)  Using  the  candidate  range  and  amplitude  estimate,  calculate  L 

5)  Repeat  Steps  2-4  until  all  candidate  ranges  have  been  tested 

6)  Select  the  range  that  corresponds  to  the  maximum  L  value 

7)  Go  back  to  Step  1  for  all  pixels  in  detector  array. 

4.3  Normalized  Cross-Correlation 

In  order  to  mitigate  inter-sample  targets,  scaling,  and  waveform  truncation  issues, 
sub-sample  ranging  is  performed  on  a  pixel’s  pulse-shape  pk  (m,  n )  (e.g.  Equation  (2.72) 
or  (2.73))  by  using  a  normalized  cross-correlation  (NCC)  method  based  on  the  Pearson 
product-moment  correlation  coefficient.  Using  this  coefficient  forces  each  pixel’s  waveform 
to  be  zero  mean  and  unit  standard  deviation.  A  symmetrical  waveform  is  assumed  for 
notation  simplicity  .  However,  an  asymmetrical  waveform  method  could  be  implemented. 

The  correlation  matrix  would  then  be  increased  by  one  dimension  due  to  breaking  up  the 
pulse-width  standard  deviation  into  two  variables:  leading  and  trailing. 

Analogous  to  a  cross-correlation  range  estimator  in  [63] ,  the  normalized  cross-correlation 
method  is  constructed  as  follows:  The  range  vector  of  samples  within  a  cube  is  represented 
by 

R  (k)  =  zmin+z[nc  (k)  (4.10) 
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where  k  E  [0, K  —  1],  K  is  total  number  of  samples  (same  K  as  defined  in  the  data 
model  in  Section  2.5),  zmm  is  the  range  of  the  first  sample,  and  zinc  is  the  range  increment 
per  sample.  Another  range  vector,  Kr  is  constructed  with  the  same  maximum  and  mini¬ 
mum  extents  as  R,  but  with  a  smaller  range  increment  per  sample  defined  by  the  following 
equation: 

Kr  (q)  =  zmiD+zf  ( q )  (4.1 1) 


where  q  E  [0,  K'  —  1],  K'  is  the  number  of  samples  in  Kr,  and  Zf  is  the  range  increment. 
Since  the  extents  of  Kr  match  R,  the  following  inequalities  hold:  K'  >  K  and  Zf  <  z-mc.  A 
2D  reference  Gaussian  waveform  matrix  is  used  with  the  Kr  vector  as  the  reference  target 
location  or 


rk  (q)  =  exp 


-(4 


2A;  jq)  /c)2 

2^ 


(4.12) 


where  4  =  2 R  (/::)  / c  and  is  the  time  vector,  R  (k)  is  the  range  vector  from  Equation  (4. 10), 
c  is  the  speed  of  light  in  vacuum,  and  crw  is  the  transmitted  pulse  standard  deviation.  The 
zero  mean  and  unit  variance  version  of  rk  for  all  k  E  [1, ...,  K\  and  q  E  [1, ...,  K ']  is 


on  ^  rk(q)~rk(q) 

S2  (fc,  q)  =  - -  (4-13) 

O’ r  (q) 

where  crr2  and  fk  are  the  variance  and  average  of  rk  in  the  time  dimension.  Considering 
the  range  estimate  for  the  (m,n)th  pixel,  the  zero  mean  and  unit  variance  version  of  the 
pulse-shape  of  interest  pk  (m,  n)  for  all  k  E  [1, ...,  K]  is 


Si  (k) 


Pk  (m,  n) 


K 

Y,  Pk(m>n) 
k=  1 _ 
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a2  (m,  n) 


(4.14) 


where  cr2  is  the  variance  of  pk  (m,  n)  in  the  time  dimension  respectively.  With  S)  and  S2 
determined,  the  normalized  cross  correlation  denoted  by  *  is  performed  by 
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The  cross  correlation  *  operation  is  carried  out  using  a  matrix  multiply  given  by 

S2  *  Si  =  (S2)T  x  S'!  (4.16) 

where  S2  and  .S')  have  dimensions  [ K ,  K ']  and  [K]  respectively,  “T”  is  the  transpose  op¬ 
erator,  and  x  is  a  matrix  multiply.  The  result  of  the  matrix  multiply  is  a  vector  of  values 
Ckt  with  dimension  [K'\  that  correspond  to  the  strength  of  the  similarity  between  the  ref¬ 
erence  waveform  S2  at  different  target  ranges  and  the  data  waveform  Si.  Finding  the  range 
estimate  is  accomplished  by  peak  detection  (i.e.  selecting  the  target  range  with  the  highest 
value  from  the  matrix  multiply)  on  C'xr  or 

R  (m,  n)  =  argmax  CKr  (q) .  (4.17) 

Zmin  +Zf  (?) 

The  NCC  method  is  used  exclusively  in  Chapter  V. 

4.4  Two  Point  Target  Range  and  Spatial  Separation  Estimator 

With  FOliage  PENetration  (FOPEN)  applications,  this  section  develops  a  range  sep¬ 
aration  estimator  by  using  a  least  squares  approach  adapted  from  previous  work  that  only 
considered  two  targets  within  a  single  pixel  in  a  non-blurry  environment  [5].  While  no 
noise  source  is  specified  in  the  subsequent  development,  estimator  results  in  a  shot-noise 
limited  environment  are  given  in  Section  6.2. 

4.4.1  Two  Point  Target  Data  Model.  The  mean  of  the  observations  in  units  of 
photons  of  a  two  point  target  scene  interrogated  by  a  3D  FLASH  LADAR  are  defined  by 
a  convolution  between  the  object  and  the  system  point-spread-function  (PSF)  added  to  a 
pixel  bias  or  [25],  [37] 

M  N 

ik  (x,  y)  =  ^2  ^2  °k  (m’ n)  -  m:V  -  n)  +  B  (x,  y )  (4.18) 

777.-1  77.— 1 
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where  (x,  y )  are  the  pixel  plane  coordinates  with  x  E  [1,  X]  and  y  E  [1,  Y],  k  is  the  range 
dimension  coordinate,  and  (m,n)  are  the  object  plane  coordinates  with  m  E  [1 ,  M]  and 
n  E  [1,  N],  The  integer  range  dimension  variable  k  E  |0.  A'  —  1]  corresponds  to  a  range 
distance  rk  in  units  of  meters  according  to 


rk  =  Ks+(^^j  (4.19) 

with  Ks  being  the  initial/starting  range  of  the  data  cube,  ts  as  the  range  sampling  interval 
in  seconds,  and  c  being  the  speed  of  light  in  meters  per  second. 


Figure  4.1:  (a)  For  illustrative  purposes,  this  figure  is  a  range  image  of  the  truth  data 

where  the  reference  target  is  in  the  center  of  the  array  at  1000  meters  with 
the  unknown  target  placed  at  Am  =  2  pixels  and  Ak  =  1.7  meters. 

(b)  Defined  by  Equation  (4.22),  this  shows  the  ideal  waveforms  of  the  un¬ 
known  p  (rk  —  Kt)  and  reference  target  p  {rk  —  Kr )  from  Figure  4.1(a)  with 
a  pulse-width  standard  deviation  apt  =  0.88  ns. 


Considering  both  range  and  spatial  dimensions,  the  two  point  target  scene  consists 
of  one  target  at  a  known  position  and  one  target  at  an  unknown  position.  The  targets  are 
constructed  this  way  since  the  paper’s  focus  is  on  range  separation  between  the  targets  and 
not  absolute  range.  This  assumption  keeps  the  parameter  of  interest  (range  separation)  in¬ 
tact  while  simplifying  the  data  model  by  preventing  an  additional  unknown  parameter.  The 
targets  are  considered  point  targets  spatially,  but  do  provide  a  returned  waveform.  Consid¬ 
ering  the  two  point  target  scene  illustrated  by  Figures  4.1(a)  and  (b) ,  the  object  ok  (x,  y)  is 
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defined  by 


ok  (m,  n)  =  Atp  (rk  -  (Kr  -  Ak))  6  (m  -  Am,  n)  +  Arp  (rk  -  Kr)  5  (m,  n) .  (4.20) 

where  At  and  Ar  are  the  point  target  amplitudes,  p  (rk  —  (Kr  —  Ak))  and  p  {rk  —  Kr)  are 
the  received  pulse  shapes  with  Kr  as  the  known  reference  target  and  Ak  as  the  range  sep¬ 
aration  between  the  known  and  unknown  target  (Kt)  or  Ak.  =  Kr  —  Kt.  While  the  range 
sampling  capability  rk  of  the  LADAR  is  fixed  by  the  receiver  electronics,  the  unknown 
target  Kt  could  occur  anywhere  within  the  range  gate  to  include  ranges  between  samples. 
Also,  the  spatial  point  targets  are  defined  by  Kronecker  delta  functions  S  (m  —  Am.  n)  and 
5  (m,  n )  and  h  (x,  y )  is  the  known  system  PSF.  The  final  term  is  the  pixel  bias  B  ( x ,  y) 
and  is  intended  to  account  for  any  ambient  light,  dark  current,  electron  noise,  and  pixel-to- 
pixel  impulse  response  variations.  This  bias  is  assumed  known  and  to  be  governed  by  the 
Poisson  distribution  due  to  the  discrete,  random  nature  of  these  noise  sources.  Concerning 
the  validity  of  the  assuming  a  known  pixel  bias,  it  is  target  independent  and  can  be  sepa¬ 
rately  determined  during  LADAR  operation  by  a  calibration  step  where  the  data  is  collected 
without  activating  the  laser. 

Performing  the  convolution  in  Equation  (4.18)  results  in  the  simplified  form 


ik  (x,  y )  =  Atp  (rk  -  (Kr  -  Ak))  h(x-  Am,  y)  +  Arp  (rk  -  Kr)  h  (x,  y)  +  B  (x,  y) 


(4.21) 


where  the  received  pulse  shapes  are  assumed  symmetric  Gaussian  and  defined  by 


(4.22) 


with  apd  as  the  pulse-width  standard  deviation  in  units  of  meters  and  defined  as  apd  = 
cxjpt/ 2  where  apt  is  the  pulse-width  standard  deviation  in  units  of  seconds.  Gaussian- 
shaped  pulses  are  a  valid  approximation  for  the  pulse  shapes  transmitted  from  3D  FLASH 
LADAR  hardware  [39].  After  analysis  on  experimental  data,  it  was  also  found  that  the 
received  pulse-shapes  display  an  inherent  asymmetry.  In  other  words,  the  pulse-shape 
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definition  is  changed  such  that  there  are  two  pulse-shape  standard  deviations  concerning 
Equation  (4.22):  one  for  the  leading  edge  (pre-target)  and  another  for  the  trailing  edge 
(post-target).  Although  the  effects  of  asymmetrical  pulses  on  the  CRB  and  range  sepa¬ 
ration  estimation  is  a  source  of  additional  research,  the  symmetry  or  lack  thereof  in  the 
received  pulses  does  not  change  the  conclusion  that  an  optimal  pulse  exists  given  the  range 
resolution  metric.  Symmetrical  pulse-shapes  are  assumed  for  simplicity  and  are  simply  a 
subset  of  asymmetrical  pulse-shapes. 

Furthermore,  a  spatially,  invariant  2D  Gaussian  PSF  is  chosen  because  its  differenti¬ 
ation  is  straight-forward  while  still  providing  a  function  to  sufficiently  blur  a  target  scene. 
This  type  of  impulse  response  has  been  used  previously  to  describe  blurring  due  to  atmo¬ 
spheric  turbulence  [37].  The  PSF  is  defined  as 


(4.23) 


where  ah  >  0  is  the  PSF  standard  deviation  (measured  in  units  of  pixels)  and  is  affected  by 
light  diffraction  effects,  receiver  optic’s  quality,  and  atmospheric  turbulence. 

4.4.2  Estimator  Derivation.  Given  the  variable  definitions  from  Equations  (4. 1 8)- 

(4.23),  the  sum  squared  error  term  is  the  sum  of  the  square  of  the  difference  between  the 
observed  data  and  the  estimate  or 


(4.24) 


k  x  y 


where  the  dependence  on  A/,,  will  subsequently  be  dropped  for  conciseness.  There  are 
four  unknowns  including  the  two  amplitudes,  range  separation,  and  spatial  separation.  The 
procedure  to  find  the  range  and  spatial  separation  estimates  is  to  iteratively  step  through 
each  possible  combination  of  range  and  spatial  separation  values  (these  values  are  known 
a  priori)  and  then  determine  the  amplitudes  that  minimize  the  error.  After  an  exhaustive 
search  of  combinations  of  range  and  spatial  separations,  the  combination  that  results  in  the 
least  sum  square  error  is  chosen  as  the  estimates. 
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For  a  particular  amplitude  Ai,  the  approach  is  to  take  the  derivative  of  the  error  (Equa¬ 
tion  (4.24))  with  respect  to  that  amplitude,  set  the  result  equal  to  zero,  dE/dAi  =  0,  and 
solve  for  the  amplitude  term.  This  method  gives  the  amplitude  value  that  minimizes  the 
error  term  due  to  the  positivity  of  the  second  derivative.  Since  there  are  two  amplitudes,  a 
well-posed  system  of  equations  is  set  up  by  performing  dE / dAt  =  0  and  dE / dAr  =  0  and 
solving  for  At  and  Ar  respectively  resulting  in  two  equations  and  two  unknowns  shown  by 

CnAt  +  Ci2Ar  =  D\ 

C2lAt  +  C*22^4r  =  D2  (4.25) 


where 

Ci i  =  (Kr-  Ak))h(x-  Am,y)}2 

k  x  y 

C\2  =  C21  =  -  (Kr  -  Ak))h(x  -  Am,y)p(k  -  Kr)h(x,y) 

k  x  y 

°22  =  Kr)h(x,y)}2 

k  x  y 

*  =  EEE  (B  ( x ,  y)  -  dk  0,  y))p  (k  -  (Kr  -  Ak))  h(x-  Am,  y) 

k  x  y 

*  =  EEE  (B  (x,  y)  -  dk  (x,  y))p  ( k  -  Kr)  h  (. x ,  y) .  (4.26) 

k  x  y 

The  amplitudes  are  then  determined  by  solving  the  system  of  equations. 

The  following  provides  the  estimation  steps: 

1 .  Select  a  range  separation,  Ak 

2.  Select  a  spatial  separation,  Am 

3.  Determine  the  estimates  for  amplitudes,  At  and  Ar,  via  the  system  of  equations 
in  (4.25) 

4.  Calculate  error  term  using  Equation  (4.24) 

5.  Repeat  steps  1-4  until  all  range  and  spatial  separations  have  been  selected 
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6.  Select  range  and  spatial  separation  corresponding  to  the  smallest  error 

4.5  Pixel-Dependent  Two  Surface  Range  Separation  Estimator 

The  purpose  of  this  estimator  is  to  be  able  to  estimate  the  range  to  two  surfaces  within 
a  single  pixel.  Similar  to  [5],  this  scenario  differs  in  that  the  first  surface  in  range  is  known 
while  the  second  surface,  further  in  range,  is  unknown.  Further,  unlike  the  previous  section 
where  the  data  model  includes  all  pixels  and  range  samples,  this  estimator  operates  on  one 
pixel  at  a  time.  Section  6.3.2  uses  this  estimator  against  a  complex,  two  surface  target  to 
find  the  optimal  pulse-width  with  respect  to  range  resolution. 

The  data  model  for  an  arbitrary  pixel  is 

dk  =  ik  +  nk  (4.27) 

where  dk  is  the  observed  data,  ik  is  the  blurry,  non-noisy  data,  and  nk  is  the  noise.  Since 
a  pixel  can  follow  more  than  one  data  model,  hypothesis  testing  is  performed  to  decide 
whether  there  is  one  or  two  surfaces  in  the  pixel.  In  a  two  surface  target  scene,  a  par¬ 
ticular  pixel  might  contain  one  or  two  surfaces.  Therefore,  the  blurry,  non-noisy  data  is 
hypothesized  to  be  either  a  “two  surface  pixel”  by 


ik  =  Atp  ( k  —  Kt )  +  Arp  ( k  —  Kr )  +  B  (4.28) 

or  “one  surface  pixel”  by 

4  =  Agp  {k  -  Kg)  +  B.  (4.29) 

where  At  and  Ar  are  the  received  target  amplitudes  (includes  convolution  effects),  p  (k) 
is  the  received  pulse  with  Kt  as  the  unknown,  second-surface  target  range  and  Kr  as  the 
known,  first-surface  target  range,  and  B  as  the  pixel  bias  term.  Note,  either  the  first  OR 
second  surface  can  represented  by  the  “one  surface  pixel”  (Equation  (4.29))  case  where  Ag 
and  Kg  are  generic  variables  representing  either  surfaces  particular  amplitude  and  range 
respectively.  The  unknown  target  range  Kt  can  also  be  defined  by  a  range  separation  Ak 
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from  the  known  target  Kr  by  Kt  =  Kr  —  Ak.  The  terms  that  must  be  estimated  (i.e.  the 
unknowns)  are  At,  Ar  (if  “two  surface  pixel”),  and  Ak. 

Regardless  of  the  number  of  surfaces  in  the  pixel,  the  sum  squared  error  metric  is  the 
same  ( E1  (Afc)  and  E2  (A/,)  for  a  “one  surface  pixel”  or  a  “two  surface  pixel”  respectively) 
and  is  defined  as 

E  (Afc)  =  (4  -  ik?  (4.30) 

k 

which  is  similar  to  the  previous  previous  section  with  the  exception  that  the  pixel  detector 
dimensions  are  dropped.  The  procedure  to  estimate  the  unknown  parameter  of  interest, 
range  separation  A^,  is  to  estimate  the  range  separation  using  both  “one  surface  pixel”  and 
“two  surface  pixel”  data  models  and  choose  the  “one  surface  pixel”  case  and  corresponding 
range  separation  estimate  if 

ZEmin 

-4—  <  7  (4.31) 

Emin  '  v  7 

where  7  is  a  threshold,  Ej""1  =  arg max  E\  (A/,.),  and  E!J"n  =  arg max  E2  (A/,.).  If  the 

Afc 

inequality  in  Equation  (4.31)  is  not  true,  then  choose  the  “two  surface  pixel”  case  and  its 
range  separation  estimate. 

The  amplitude  and  range  estimate  on  the  “one  surface”  data  model  are  attained  by 
selecting  a  candidate  range  separation  and  then  taking  dEi/dAt  =  0  and  then  solving  for 
At  given  by 

E  {dkp  (. k  -  Kt)  -  Bp  (k-  Kt)} 

At  =  — - - - .  (4.32) 

Ep(k-Kt)2 

k= 1 

The  “two  surface”  amplitude  estimates  are  determined  in  the  same  manner  as  the 
previous  sections  using  Equations  (4.25)  and  (4.26)  and  dropping  the  pixel  dimension  ( x 
and  y)  summations. 

The  steps  of  the  estimator  are: 

1 .  Select  a  range  separation,  Ak 


82 


2.  Determine  the  amplitude  estimates  for  the  “two  surface  pixel”  case  via  Equations  (4.25) 
and  (4.26) 

3.  Determine  the  amplitude  estimate  for  the  “one  surface  pixel”  case  via  Equation  (4.32). 

4.  Calculate  error  terms  E\  and  E2  using  Equation  (4.30) 

5.  Repeat  steps  1-4  until  all  range  separations  have  been  selected 

6.  Using  the  hypothesis  test  from  Equation  (4.31),  select  the  range  separation  corre¬ 
sponding  to  the  smallest  error.  The  range  separation  for  the  pixel  may  be  zero  if  the 
“one  surface  pixel”  case  is  chosen  and  if  the  pixel  is  a  “first  surface  pixel”  as  well 
with  a  known  range  of  Kr. 
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V.  Improving  3D  FLASH  LADAR  Range  Estimation  via  Object 

Recovery 


The  motivation  in  this  chapter  is  to  provide  a  means  of  improving  range  estimation  by 
object  recovery  (i.e.  spatially  deblurring  data)  from  three  dimensional  FLASH  LAser 
Detection  And  Ranging  (3D  FLASH  LADAR)  observations.  Referring  to  Figure  5.1(a), 
the  idea  is  to  process  the  data  in  the  spatial  dimensions  (x,  y )  while  improving  ranging 
performance  in  the  time  dimension  (k).  Taken  exclusively  from  [56],  this  chapter  covers 
novel  material  including  amplitude,  pulse-shape,  system  impulse  response,  and  pixel  bias 
estimation.  Original  efforts  also  include  object,  system  impulse  response,  and  pixel  bias 
estimation. 

Building  on  material  presented  in  Chapters  II  and  III,  a  method  to  model  the  3D 
FLASH  LADAR  data  operating  in  “sular  mode”  is  that  the  2D  range  images  are  formed 
via  a  convolution  between  the  object  at  a  particular  time  and  the  spatial  impulse  response. 
In  Figure  5.1(a),  a  range  image  d(tk)  is  one  of  the  2D  slices  of  the  data  cube.  Considering 
the  laser  illuminating  a  target,  one  collect  from  a  3D  FLASH  LADAR  sensor  results  in  a 
data  cube  consisting  of  a  series  of  range  images  ( N  from  Figure  5.1)  representing  detected 
photons.  (NOTE:  Figure  5.1  is  shown  again  for  convenience.) 

Attempts  at  3D  FLASH  LADAR  range  estimation  of  a  remote  scene  can  result  in 
errors  due  to  several  factors  including  the  optical  spatial  impulse  response,  detector  blur¬ 
ring,  photon  noise,  timing  jitter,  and  readout  noise.  These  factors  either  cause  the  scene’s 
intensity  to  spread  across  pixels  or  add  unwanted  and  disruptive  noise  effects.  The  intensity 
spreading  and  noise  corrupt  the  correct  pixel  intensities  by  mixing  intensities  with  neigh¬ 
boring  pixels  thereby  providing  false  intensity  values  and  therefore  incorrect  photon  counts 
to  the  range  estimator.  Without  blur  and  noise  compensation,  the  range  estimates  would  be 
inaccurate  to  a  degree  depending  on  the  blur  and  noise  severity. 

The  theoretical  development  of  the  range  estimator  algorithm  is  covered  first  and 
then  verified  using  simulation  and  experimental  results.  The  algorithm  is  a  variation  on  the 
Expectation  Maximization  (EM)  algorithm  called  Generalized  Expectation  Maximization 
(GEM)  and  is  desirable  due  to  its  iterative  likelihood  maximization,  convergence  proper- 
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Figure  5.1:  (a)  3D  view  of  LADAR  system  model  in  Cartesian  coordinates  with  each 

data  cube  having  dimensions  of  pixel  x  pixel  x  time  sample.  The  variable 
d(tk)  corresponds  to  the  kth  receiver  detected  range  image  with  k  G  [1,  N]  . 
(b)  Another  view  of  the  3D  FLASH  LADAR  operation.  Each  range  image’s 
full  field  of  view  (FOV)  is  128  x  128  pixels  with  a  range  gate  near  2  nanosec¬ 
onds  corresponding  to  the  3D  FLASH  LADAR  system  used  for  experimental 
collects. 
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ties,  and  ability  to  decouple  terms  [54].  The  GEM  algorithm  is  powerful  in  that  it  can 
perform  blind  deconvolution  in  situations  with  severe  defocus  or  other  aberrations  includ¬ 
ing  atmospheric  turbulence.  To  account  for  different  scenarios,  two  versions  of  the  GEM 
algorithm  are  derived  that  either  recover  the  pulse-shape  or  the  object.  The  primary  differ¬ 
ence  between  the  two  involves  data  required  and  accuracy.  Pulse-shape  estimation  requires 
less  data,  but  is  less  accurate  than  object  estimation.  Additional  details  of  the  differences 
are  presented  in  Sections  5.1.2  and  5.1.3.  Both  pulse-shape  and  object  GEM  algorithms  are 
novel  contributions  to  the  research  area. 

In  addition  to  the  GEM  algorithms,  a  Wiener  filter  method  is  used  to  attempt  range  es¬ 
timation  improvement  via  object  recovery  from  3D  FLASH  LADAR  observations  [37] ,  [55] . 
Requiring  spatial  impulse  response  knowledge  a  priori,  this  method  can  only  perform  de- 
convolution  unlike  the  blind  deconvolution  ability  of  the  GEM.  The  purpose  for  adding 
this  other  method  is  to  show  that  the  GEM  outperforms  a  competing  algorithm  that  already 
knows  part  of  the  answer  (spatial  impulse  response). 

This  chapter  is  organized  as  follows:  Section  5.1  describes  the  Wiener  filter  theory 
and  derives  the  pulse  and  object  iterative  estimators  via  the  GEM  algorithm,  Section  5.2 
presents  results  from  simulated  data  showing  improvement  in  range  estimation  after  object 
recovery,  Section  5.3  details  the  results  from  an  experimental  collection  and  processing 
also  showing  an  improvement  in  range  estimation  after  object  recovery,  and  Section  5.4 
provides  conclusions  based  on  observed  results. 

5.1  Theoretical  Development 

This  section  details  the  object  recovery  methods  used  on  simulated  and  experimen¬ 
tal  data.  Even  though  the  laser  light  is  partially  coherent,  the  argument  is  made  that  the 
detected  light  is  able  to  be  modeled  as  fully  incoherent.  The  incoherent  light  model  still 
captures  the  discrete,  non-negative  nature  of  the  received  photons  that  the  partially  coher¬ 
ent  model  exhibits.  In  addition,  experimental  data  processing  from  Section  3.6  showed  that 
the  speckle  parameter  estimation  results  tend  towards  the  incoherent  model.  Consequently, 
this  incoherent  light  model  assumption  allows  for  the  returns  to  be  a  result  of  a  linear, 
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spatially  invariant  (LSI)  system  involving  an  intensity  convolution  (instead  of  amplitude 
convolution)  between  the  intensity  point  spread  function  and  the  remote  scene.  Linearity  is 
a  consequence  of  electromagnetic  wave  propagation  theory,  and  spatial  invariance  results 
from  remaining  with  the  isoplanatic  angle  [25].  Utilizing  this  LSI  convolution  model,  two 
GEM  blind  deconvolution  algorithms  are  derived  that  enable  improved  range  estimation. 
All  references  to  the  scenario  or  data  model  refer  to  material  presented  in  Section  2.5. 

The  unknown  parameters  in  this  scenario  are  the  object  (target  amplitude  and  target 
range)  PSF,  and  bias.  The  variable  of  interest  in  this  paper  is  the  range  term  residing  in 
Equation  (2.72)  or  (2.73).  Direct  estimation  of  the  range  term  is  problematic  because  of 
its  location  either  in  an  exponential  or  in  a  squared  term.  Therefore,  the  approach  to  range 
estimation  is  to  retrieve  the  range  from  the  estimated  pulse-shape  or  object.  This  method¬ 
ology  relies  on  the  knowledge  that  the  target  produces  the  waveform  peak  in  the  detected 
returns.  Concerning  the  PSF,  blind  deconvolution  techniques  must  be  employed  since  the 
PSF  is  unknown.  Blind  deconvolution  has  a  rich  heritage  in  astronomical  imaging  provid¬ 
ing  a  bevy  of  literature  attempting  blind  deconvolution.  Although,  blind  deconvolution  in 
astronomical  cases  consists  of  trying  to  recover  one  object  and  one  PSF  (or  many  PSFs  if 
using  multiple  frames).  In  trying  to  recover  the  target  range  from  one  3D  FFASH  FADAR 
data  collect,  this  problem  consists  of  many  objects  with  one  PSF.  There  are  many  objects 
due  to  the  transmitted  waveform  causing  each  range  slice  to  contain  different  intensities 
corresponding  to  where  the  waveform  is  incident  on  the  object.  Therefore,  these  incident 
points  become  distinct  objects  in  the  blind  deconvolution  framework.  If  multiple  cubes 
are  necessary,  the  atmosphere  is  changing  with  each  cube  resulting  in  multiple  PSFs  that 
must  be  estimated  resulting  in  a  “multi-frame”  or  “multi-cube”  scenario.  If  no  atmospheric 
turbulence  exists  or  is  non-volatile,  the  PSF  is  consistent  throughout  the  cubes  and  the  j 
subscript  can  be  dropped. 

5.1.1  Object  Deconvolution.  As  noted  previously,  the  goal  of  this  research  effort 
is  to  improve  range  estimation  of  a  target  illuminated  by  a  3D  FFASH  FADAR.  Deconvo¬ 
lution  is  necessary  due  to  the  imaging  nature  of  the  3D  FADAR  producing  blurred  return 
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pulses.  A  solution  is  to  use  image  restoration  techniques  to  attempt  to  restore  the  original 
range  slice  images  thereby  improving  the  range  estimation.  The  image  restoration  algo¬ 
rithm  applies  a  2-D  filter  to  the  a  pixel  detector  array  at  each  time  sample  (i.e.  range  slice 
image,  etc.)  resulting  in  a  “de-blurred”  data  cube.  The  “de -blurred”  data  cube’s  pixels  now 
more  closely  mimic  the  unblurred  return  pulse  from  Equation  (2.72)  and  result  in  improved 
range  estimation. 

A  standard  linear  filter  that  can  perform  the  image  restoration  is  the  pseudo-inverse 
Wiener  filter.  From  [37],  the  definition  of  the  pseudo-inverse  Wiener  filter,  Gw(fx,  fy)> 's 


Gw(fx j  fy )  — 


\H(fXJy)\*  + 


1 


SNR 


(5.1) 


where  H(fx ,  fy)  is  the  overall  optical  transfer  function  (OTF),  *  is  the  conjugate  operator, 
and  SNR  is  the  signal-to-noise  ratio  in  the  image.  One  image  processing  definition  of  the 
SNR  is  to  set  it  equal  to  the  image  mean  divided  by  the  image  standard  deviation  [70]. 
Given  that  the  signal  is  dominated  by  shot  noise,  the  SNR  is  defined  at  particular  range 
sample  k  to  be  the  mean  of  the  detected  range  image  pd  divided  by  the  detected  range  image 
standard  deviation  ^/Jif  or 

SNR  =  / JTd .  (5.2) 

Using  the  pseudo-inverse  Wiener  filter,  the  deblurred  image  at  a  particular  range  sample  k 
is 

h{x,  y )  =  F~l  {Gw  (fx,  fy)  Dk  (fx,  fy)}  (5.3) 


with  F  1  as  the  inverse  Fourier  transform  and  Dk  (fxi  fy)  is  the  Fourier  transform  of  the 
detected  range  image  dk(x,y).  After  the  cube  has  been  filtered,  the  normalized  cross¬ 
correlation  (NCC)  range  estimator  method  uses  the  filtered  data  to  determine  the  range 
estimate  using  Equation  (4.17).  The  waveform  variable  pk  (m,n)  in  Equation  (4.14)  is 
replaced  by  Ik(x,  y)  during  the  NCC  implementation. 


5.1.2  Pulse-Shape  Blind  Deconvolution  via  the  GEM  Algorithm.  The  previous 
section  assumed  a  known  PSF.  This  section  covers  the  blind  deconvolution  case  where  the 


88 


PSF  is  unknown.  Blind  object  recovery  is  accomplished  using  two  approaches  concerning 
the  pulse-shape  and  object  variables  from  Equation  (2.72).  The  pulse-shape  estimation  is 
very  powerful  in  that  the  estimator  only  needs  one  data  cube  (one-shot,  one-kill).  However, 
if  the  best  accuracy  is  required  and  the  3D  data  cubes  are  properly  registered,  the  multi-cube 
object  estimation  provides  lower  error. 

Referring  to  the  GEM  theory  from  Section  2.4  and  the  data  model  from  Section  2.5, 
the  model  is  reformed  to  consider  pulse-shape  recovery  with  one  cube  required  for  process¬ 
ing  (with  j  =  1).  Consistent  with  the  GEM  algorithm,  the  original  data  dk(x,  y )  is  called 
the  incomplete  data  and  is  defined  by  [54] 


M  N 


(5.4) 


777—1  77=1 


where  two  new  variables,  dk(x,  y\m,  n )  and  qk  (x,  y ),  are  called  complete  data.  This  formu¬ 
lation  provides  two  sets  of  complete  data  that  account  for  the  photon  noise/image  formation 
and  pixel  bias  respectively.  The  formation  of  the  complete  data  highlights  the  powerful  na¬ 
ture  of  the  GEM  algorithm.  In  this  application,  complete  data  can  also  be  called  unobserved 
data  and  carries  no  explicit  physical  meaning.  It  is  used  to  directly  benefit  the  EM  algo¬ 
rithm.  Consistent  with  [71],  careful  definition  of  the  complete  data  allows  decoupling  of 
unknown  variables  while  preserving  physical  meaning  in  the  expected  value  of  the  incom¬ 
plete  data. 

The  expected  value  of  the  complete  data  sets  is  given  by 


E  dk(x,y\m,n)  =  A(m,n)  pk  (m,n)  h  (x  —  m,y  —  n) . 


(5.5) 


and 


E  [qk  (x,  y)}  =  B  (x,  y) 


(5.6) 


where  B  (x,  y )  is  the  constant  pixel  bias.  The  expected  value  of  the  incomplete  data  is  thus 


E  [dk  (x,  y)]  =  ik  (x,  y)  +  B  (x,  y) 


(5.7) 
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which  is  consistent  with  the  data  observations  depicted  in  the  high  fidelity  model  (i.e.  con¬ 
volution  model)  of  Figure  2.4.  Adding  the  pixel  bias  to  the  model  covers  non-modeled 
noise  effects  and  pixel-to-pixel  impulse  response  variations.  The  pixel  bias  is  assumed  to 
be  governed  by  the  Poisson  distribution  due  to  the  discrete  random  nature  of  dark  current 
and  electron  noise.  Physically,  the  pixel  bias  is  added  to  the  photons  incident  upon  the 
detector  and  is  part  of  the  detected  photon  counts.  The  PMF  for  the  photon  noise  is 


P  (dk  (x,y\m,n)^j  = 

[A  (m,  n )  pk  (m,  n)  h  (x  —  m,y  —  n)]^k^x,y^  Q-iA(m’n)Pk(m,n)h(x-m,y-n)} 


dk  {x,y)\ 


(5.8) 


while  the  PMF  for  the  pixel  bias  is 


P(qk  (x,y))  = 


B(x,y)gAx,y^e  B(X'A 

Qk(x,y)\ 


(5.9) 


Assuming  statistical  independence  between  all  the  pixels  and  between  the  photon  noise  and 
pixel  bias  noise,  the  complete  data  log-likelihood  function  considering  all  random  variables 
is 


Lcd  ( Pk ,  4  h,  B )  =  In 


n 


P  (dk  (x,y\m,n))P  (qk  (x,y)) 


_k,x,y,m,n 

or  (NOTE:  summations  wrap  around  unless  otherwise  stated) 


(5.10) 


Lcd  (Pk,  4  h,  B)  =  ^  dk  (x,  y\m,  n)  In  [A  (m,  n)pk  (m,  n)  h  (x  -  m,  y  -  n)] 

k,x,y,m,n 

-  [A  (m,  n)  pk  (m,  n)  h  (x  -  m,  y  -  n)}  +  qk  (x,  y)  In  [B  (x,  y)\  -  B  (x,  y) .  (5.11) 


Referring  to  Equation  (2.64),  the  Q  function  then  becomes 


Q  ( Pk ,  4  h,B)  =  E  [Lcd  ( Pk ,  A,  h,  B)  |4  (x,  y) , pf ,  A°ld,  h°ld ,  B°ld ] 


(5.12) 
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where  the  estimates  for  the  amplitude,  pulse-shape,  PSF,  and  bias  are  considered  maximum- 
likelihood  estimates.  Taking  the  conditional  expectation  of  Equation  (5.12)  results  in 


Q  ( Pk ,  A,  h,  B )  = 

(x’  y> m’ n;  A<*d' pkd’  h°ld )  ln  \-A  (m’ n) Pk  (m’  n)h(x-m,y-  n)} 

k,x,y,m,n 

-  [A  (m,  n )  pk  (m,  n )  h(x-m,y  -  n)}  +  y°±d  (x,  y\ Bold )  In  [5  (x,  y)] 
-B(x,y)  (5.13) 


where 


(xiyimin\P0kd,Aold,hold)  =  E  \dk  (x,y\m,n)  \dk(x,y) , p°kld ,  Aold ,  h 


, old 


(5.14) 


and 


(i,  V,  B°M)  =  E  [qk  (x,  y)  \dk  (x,  y) ,  B°ld]  . 


(5.15) 


Equations  (5.14)  and  (5.15)  represent  the  expected  value  of  one  set  of  complete  data  given 
the  incomplete  data.  For  Poisson  random  variables,  these  expectations  turn  out  to  be  ratios 
of  the  data  times  one  expected  value  of  the  complete  data  divided  by  the  all  sets  of  expected 
values  of  the  complete  data  [75].  For  the  first  set  of  complete  data,  dk  (x,y),  the  conditional 
expectation  is 


dk  (x,  y)  A°ld  (m,  n)  pkd  (m,  n)  hold  (x  —  m,y  —  n) 


ikd  (x,  y)  +  Bold  (x,  y) 


(5.16) 


while  the  second  set  of  complete  data  concerning  the  pixel  bias  qk  (x,  y),  has  a  conditional 
expectation  equal  to 


/4d  (*,  y\ B°ld ) 


dk  (x,  y)  Bold  (x,  y) 
i°kld  (x,  y)  +  Bold  (x,  y) 


(5.17) 


The  maximization  of  the  Q  function  for  all  parameter  unknowns  (target  amplitude,  target 
pulse  shape,  PSF,  and  pixel  bias)  is  theoretically  intractable  due  to  coupling.  It  is  required 
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to  break  apart  the  Q  function  into  separate  components  such  that 


Q  —  Qp  +  Qh  +  Qa  +  Qb 


(5.18) 


where  each  component  of  the  Q  function  can  be  maximized  independently.  Thus,  the  GEM 
algorithm  becomes 


QP{pnkew\Pkd,A°ld,h°ld) 
Qa  ( Anew\p°kld,Aold ,  hold ) 
Qh  ( hnew\p°kld,Aold ,  hold ) 
Qb  (. Bnew\Bold ) 


>  gft  (hold\p°kld,Aold1hold) 

>  QB(Bold\Bold ) 


(5.19) 


which,  if  these  conditions  are  met,  ensures  that  the  likelihood  is  increased  with  each  itera¬ 
tion  [54] 

L  (p™w,  Anew,  hnew ,  Bnew)  >  L  (p°kld,  Aold ,  hold ,  5oW)  (5.20) 

resulting  in  a  GEM  sequence  converging  to  a  local  maximum. 

Beginning  the  estimation  process  of  the  separate  Q  functions  starts  with  the  target 
pulse  shape,  Qp  which  is 

'  K 

Qp=  X/  ^4  Og  20  w;  Pfc  d,  hold)  In  [pfe  (m,  n)]  -  A  (m,  n)  (m,  n)  -  1 

k,x,y,m,n  \_  k=l 

(5.21) 

where  a  pixel-dependent  Lagrange  multiplier,  A  (m,n),  is  introduced  to  force  the  pulse 
shape  to  add  to  one  for  each  pixel.  This  constraint  is  necessary  to  decouple  the  pulse 
shape  from  the  target  amplitude  and  PSF.  Taking  the  derivative  of  Equation  (5.21)  with 
respect  to  a  particular  object  plane  point  and  range  sample,  setting  the  result  equal  to  zero, 
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dQp/dpk0  (m0,  n0)  =  0,  and  solving  for  the  pulse  shape,  results  in 


PkT"  (mc>)  n0)  =p°ko  ( m0,n0 ) 


Aold  (m0,  na)  \  ^  ^  40  (xi  v)  hold  (x  -  m0,y  -  n0) 


A  (m0,  n0)  )  ^  ^  (x,  y)  +  B°ld  (x,  y) 


(5.22) 


where 


A  K  ^  (m0,  nj  x:  (rn„,  n„)  E  E  4  ^  ^  ^ 

fc  =  l  X=1  J/= 


^  iokld  0,  y )  +  B°ld  ( x ,  ?/) 


(5.23) 


and 


M  N 


i°kld  (x,  V)  =  Y1  AM  ( m ’ n)  Pko  (m> n)  h°M  ( x-m,y-n ).  (5.24) 


ra=l  n=l 


Equation  (5.22)  is  the  iterative  solution  for  the  pulse  shape  per  range  sample.  Next,  the  Q 
function  is  partitioned  into  its  target  amplitude  components 


M  N 

Qa=  {pi?  {x’y’m’n'iPkld,Aold,hold)  ln[Al(m,n)]}  (5.25) 

k,x,y,m,n  m=  1  n= 1 

where 


x  Y 

^2^h(x,y)  =  l  (5.26) 

x=l  y= 1 
K 

^pk(m,n)  =  1  (5.27) 

fc=i 


have  been  utilized  to  decouple  the  pulse  shape  and  PSF  terms  from  the  target  amplitude. 
Maximizing  Equation  (5.25)  by  8Qa/9A  (m0,  n0)  =  0  and  solving  for  the  amplitude  term 
results  in  the  iterative  solution  for  the  target  amplitude  term 


K 


Ar‘ 


\jrio-,  no, 


=  A 


old 


'jn0,n0)  J2Pkd 


m0,  n0 


k= 1 


x  Y 

EE 

x=l  y=  1 


dk  (x,  y)  hold  (x  —  m0,  y  —  n0) 


i°kld  (x,  y)  +  Bold  (x,  y) 


(5.28) 
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The  PSF  is  the  final  unknown  parameter  that  uses  the  first  set  of  complete  data,  c4  (x,  y ). 
The  Q  function  for  the  PSF  is 

Qh  = 

X]  (z,  y,  m,  n;  pf,  A°ld ,  h°ld)  In  [h(x-m,y-  n)} 

k,x,y,m,n 

—  [A  (m,  n )  pk  (m,  n)  h  (x  —  m,  y  —  n)\ ,  (5.29) 

which  still  has  the  target  amplitude  and  pulse  shape  terms.  Similar  to  [71],  a  change  of 
variables  is  required  to  remove  the  dependence  on  the  pulse  shape  and  to  allow  for  easier 
differentiation.  Utilizing  J2k=iPk  (m>  n)  =  1  and  setting  m!  =  x  —  m  and  n'  =  y  —  n,  Qh 
then  becomes 


Qh  — 

{/i|“  (x,  y,x-m',y-  n';  pf,  A°ld ,  In  [/i  (m',  n')] } 

k,x,y,m'  ,n' 

—  A(x  —  m!,y  —  n')  h(mf,n').  (5.30) 

x,y,m'  ,n' 

Setting  dQh/dh  (m'0,  n'0)  =  0  and  solving  for  the  PSF  produces  the  iterative  solution 

hnew  /  /  >  S  _  hoid  /  /  /  \  V"  4  (s,  I/)  2loM  (x  -  m!0,  y  -  <)  pgM  (x  -  m'0,  y  -  n'0) 

n  rijo)  v'^o?  v  Y  y 

(x,  y)  +  5oW  (x,  y))  S  Anew  (x  —  m'0,  y  —  n'0) 

x—\  y=  1 

(5.31) 

Usually,  the  target  amplitude  term  in  the  denominator  would  be  an  issue  because  it  is  con¬ 
sidered  the  new  estimate.  However,  Equation  (5.28)  is  the  new  estimate  and  can  replace 
the  target  amplitude  in  the  denominator  resulting  in  a  consistent  solution  for  the  PSF.  Fi¬ 
nally,  the  pixel  bias  must  be  estimated.  In  order  to  estimate  the  pixel  bias,  the  second  set  of 
complete  data,  (x,  y),  is  utilized.  The  Q  function  for  the  pixel  bias  is 


K  X  Y 

Qb  =  X^X^X^7 

k= 1  x=l  y=  1 


dk  (x,  y)  Bold  (x,  y) 
i°kd  {xi  V )  +  Bold  (x,  y) 


In  (B  (x,y))  -B(x,y). 


(5.32) 
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Setting  dQs/dB  (xa,  y0)  =  0  and  solving  for  the  pixel  bias  results  in  an  iterative  solution 


(5.33) 


corresponding  to  the  pixel  bias  iteration. 

After  a  pre-determined  number  of  iterations  on  Equations  (5.22),  (5.28),  (5.31),  and 
(5.33),  range  estimate  updates  for  each  pixel  are  generated  by  using  the  NCC  method  be¬ 
tween  a  reference  waveform  at  sub- sample  ranges  and  the  the  GEM  estimate  for  the  pulse 
shape,  p^ew.  The  range-dependent  reference  waveform  that  results  in  the  highest  correla¬ 
tion  is  chosen  and  the  corresponding  range  is  the  new  range  estimate  for  that  pixel.  The 
new  range  estimate  is  fed  back  into  the  pulse-shape  to  generate  a  new  p°^d  followed  by 
another  set  of  GEM  iterations.  The  process  (GEM  iterations  followed  by  range  updates) 
repeats  with  the  new  range  estimates  used  in  calculating  p°^d  using  Equation  (5.22)  and 
ceases  when  the  mean  square  error  (MSE)  between  the  data  and  non-noisy  range  images 
reaches  the  stopping  criteria.  All  previous  amplitude,  PSF,  and  pixel  bias  estimates  carry 
over  from  one  range  update  to  the  next.  More  specifically,  iterations  cease  when  the  MSE 
is  lower  than  the  average  data  variance  or 

I<  X  Y  I<  X  Y 


Y  Y  Y  ( d » (*.  y)  -  4“‘  (*,  y)  -  Bnew  (x,  y))1  <  y  Y  Y 14  y)  (5-34) 


k= 1  x=l  y=  1 


k= 1  x=l  y=  1 


with 


M  N 


771=1  71=1 


and 


V 


(5.36) 
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where  Vk  is  the  variance  for  the  volume  elements  (voxels).  In  the  experimental  data,  a  spe¬ 
cific  distribution  for  the  variance  is  not  chosen  in  order  to  account  for  all  noise  sources.  For 
the  simulation  data,  the  primary  noise  source  is  defined  explicitly  by  the  Poisson  distribu¬ 
tion.  Therefore,  data  variance  14  is  defined  by  the  variance  of  the  Poisson  distribution  (i.e. 
mean  of  the  data). 

In  summary,  the  pulse-shape  estimation  algorithm  steps  are: 

1 .  Initialize  PSF,  amplitude,  and  pixel  bias 

2.  Determine  initial  ranges  and  define  pulse-shape 

3.  Perform  GEM  iterations  using  Equations  (5.22),  (5.28),  (5.31),  and  (5.33) 

4.  Use  NCC  to  find  new  range  estimates  with  Equation  (4.17) 

5.  Generate  new  pulse-shapes  based  on  new  ranges 

6.  Determine  MSE  and  compare  to  stopping  criteria 

7.  Repeat  Steps  3  through  6  until  stopping  criteria  violated 

8.  Range  estimates  taken  from  last  execution  of  Step  4 

In  step  one,  the  PSF  is  initialized  by  the  diffraction-limited  PSF  of  the  system  with  some 
defocus  to  allow  the  iterations  the  freedom  to  modify  the  estimate.  The  amplitudes  and  the 
pixel  bias  are  initialized  by  a  matrix  of  ones. 

5.1.3  Object  Blind  Deconvolution  via  the  GEM  Algorithm.  When  multiple  cubes 
are  available  and  properly  registered  spatially  and  temporally,  another  method  to  perform 
range  estimation  is  to  relax  the  constraint  on  the  pulse-shape  and  assume  just  an  object 
in  the  data  model.  This  change  mitigates  the  issue  in  the  hardware  data  where  the  pulse- 
shape  is  vaguely  known.  Therefore,  estimation  is  performed  on  ok  rather  than  on  pk  from 
Equation  (2.71).  The  problem  setup  is  similar  to  the  pulse-shape  estimation  (now  with  more 
than  one  cube)  by  calling  the  original  data,  djk(x,  y),  the  incomplete  data  and  specifying 


M  N 


(5.37) 


777—1  77—  1 
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where  two  new  variables,  djk(x,  y\m,  n)  and  qjk(x,y),  are  defined  and  called  complete 
data.  This  formulation  provides  two  sets  of  complete  data  that  account  for  the  image  for¬ 
mation  and  pixel  bias  respectively.  The  same  PSF  can  be  assumed  for  adjacent  collections 
due  to  a  typical  data  collection  scenario  where  environments  shouldn’t  be  changing  rapidly 
(ignore  j).  Thus,  the  expected  values  of  the  complete  data  sets  are  given  by 


E 


djk  ( x,y\m,n ) 


ok  (m,  n)  h  (x  —  m,y  —  n) 


(5.38) 


and 

E  [ qjk  ( x ,  y)]  =  B  (x,  y)  (5.39) 

where  B  ( x ,  y)  is  the  constant  pixel  bias.  The  expected  value  of  the  incomplete  data  is  thus 


E  [djk  (x,  y)}  =  ik  (x,  y)  +  B  (x,  y) . 


(5.40) 


The  PMF  for  the  photon  noise  is 


P  (dk  (x,  y\m,  n)  j  = 

[ok  ( m ,  n)  h  (x  —  m,y  —  n^dikP’y\m’n'>  eXp  {—ok  (m,  n)  h  (x  —  m,  y  —  n)} 


djk  ( x,y\m,n ) 


(5.41) 


while  the  pixel  bias  PMF  is 


P(qjk  (x,y))  = 


B(x,y)qjk^x,y^e  B(x’y ) 
Qjk  {x,y)\ 


(5.42) 


Assuming  statistical  independence  between  all  the  pixels  and  between  the  photon  noise  and 
pixel  bias  noise,  the  complete  data  log-likelihood  is  then 


Lcd  (ok,  h,  B)  —  In 


n  p[d  jk  (x,  y\m,  n)jP  (qjk  (x,y)) 

j,k,x,y,m,n 


(5.43) 
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or 


LCd  (ok,  h,  B)  =  ^2  djk  (x,y\m,n)  In  [ok  (m,  n)  h(x  -  m,y  -  n)\ 

j,k,x,y,m,n 

-  [ok  (m,  n)  h(x  -  m,y  -  n)\  +  qjk  (. x ,  y )  In  [B  (x,  y)}  -  B  (x,  y) .  (5.44) 

Referring  to  [54],  the  Q  function  becomes  the  expected  value  of  the  complete  data  log- 
likelihood  function  with  respect  to  the  incomplete  data  and  old  parameter  estimates 


Q  (Ok,  h,B)  —  E  [Lcd  (Ok,  h,  B)  \djk  (x,  y) ,  of,  hold,  Bold]  .  (5.45) 


Taking  the  conditional  expectation  from  Equation  (5.45)  results  in 

Q(ok,h,B)  =  Y'  nf  (x,  y,  m,  n;  of,  hold)  In  [ok  (m,  n)  h(x  —  m,y  —  n)} 

j,k,x,y,m,n 

-  [ofc  (m,  n)h(x  -  m,y  -  n)}  +  [if  (x,  y;  Bold )  In  [B  ( x ,  y)\  -  B  (x,  y)  (5.46) 

where 

nt l(x,y,m,n;o?,hM)  = 


and 


Similar  to  the  pulse-shape  estimation,  the  maximization  of  the  Q  function  for  all  parameter 
unknowns  (object,  PSF,  and  pixel  bias)  is  theoretically  intractable  due  to  coupling.  It  is 


E 


djk  (x,  y\m,  n)  \djk  (x,  y) ,  of,  hold 


djk  (x,  y)  of  (m,  n)  hold  {x  —  m,y  —  n) 
if  ( x ,  y)  +  Bold  ( x ,  y) 


(5.47) 


E  [qjk  (x,y)  |  djk  (. x,y),Bold ] 

djk  {x,  y)  Bold  {x,  y) 
if  (x,  y )  +  Bold  (x,  y) 


(5.48) 
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required  to  break  apart  the  Q  function  into  separate  components  such  that 


Q  —  Qo  +  Qh  +  Qb  (5.49) 

where  each  component  of  the  Q  function  can  be  maximized  independently.  Thus,  the  GEM 
algorithm  becomes 


Go  (or” 

> 

Go  (of  I 

of,  hold ) 

r\  (■knew\„old  u°ld\ 

hh  [h  \ok  ,  n  J 

> 

Qh  (ft°“ 

|  of,  hold ) 

Qb  (Bnew\Bold) 

> 

Qb  (Bold\Bt*d) 

ensuring  that  the  likelihood  is  increased  with  each  iteration  [54] 

L  {onkew,  hnew ,  Bnew)  >  L  ( o°kld ,  hold,  Bold )  (5.51) 


resulting  in  a  GEM  sequence  converging  to  a  local  minimum.  The  procedure  to  find  the 
maxima  of  the  Q  functions  is  the  same  as  in  pulse-shape  estimation.  First,  the  object  solu¬ 
tion  is  found  by  specifying 


Qo  =  ff  (x,  y,  m,  n\  o°kld,  hold )  In  [ok  (m,  n)\  -  ok  (m,  n)  h  (x  -  m,y  -  n) . 

j,k,x,y,m,n 

(5.52) 

In  order  to  maximize  Q0,  the  derivative  of  Equation  (5.52)  with  respect  to  a  particular  object 
plane  point  and  range  sample  is  set  equal  to  zero,  dQp/doko  (m0,  n0)  =  0.  Solving  for  the 
object  results  in 


°ZW  (m°’  no )  = 


^ old 
}k0 


m0:  n, 


J 


J  X  Y 

ZEE 

j=  1  X=1  y=  1 


djk0  (x,  y)  hold  ( x  —  m0,  y  —  n0) 
ikd  (x,  y)  +  B°ld  (x,  y) 


(5.53) 


99 


with  J  as  the  number  of  cubes  and  utilizing 


x  Y 

^2  '22  h  (x,  y)  =  1  (5.54) 

x=l  y= 1 

and  where 

M  N 

if  ( x ,  y)  =  2222  °kd  (m’ n)  h°ld  ( x-m,y-n ).  (5.55) 

m= 1  n=l 

Equation  (5.53)  is  the  iterative  solution  for  the  object  per  range  sample.  The  PSF  is  the 
other  unknown  parameter  that  uses  the  first  set  of  complete  data,  djk  (x.  y ).  The  0  function 
for  the  PSF  is 

Qh  =  ^2  y°Jdk  {x’y’mini  °kdi  h°ld) ln  [h(x  -  m,y  -  n)\—Ok  (m, n)  h  (x  -  m,y  -  n). 

j,k,x,y,m,n 

(5.56) 

Similar  to  [71],  a  change  of  variables  is  required  to  remove  the  dependence  on  the  pulse 
shape  and  to  allow  for  easier  differentiation.  Setting  m!  =  x  —  m  and  n’  =  y  —  n,  Qh  then 
becomes 


Qh  — 

{>JZ  y^x~m'^y~  n'’  °kdi  h°ld) ln  lh  n')] } 

j  ,k  ,x  ,y  ,m' ,n' 

— Ok  {x  —  m1  ,y  —  n1)  h  (m1  ,n .  (5.57) 


Setting  dQh/dh  (mf0,  n'0)  =  0  and  solving  for  the  PSF  produces  the  solution 


fcBewK,<)  = 


hold  (m'0,  n'0) 


J 
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(5.58) 


The  object  term  in  the  denominator  is  the  new  estimate  from  Equation  (5.53).  Since,  there 
are  phase  abberations  across  the  aperture  and  the  PSF  needs  to  be  constrained,  phase  re¬ 
trieval  is  performed  on  Equation  (5.58)  by  the  Gerchberg-Saxton  algorithm  [23].  In  the 
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pulse-shape  estimation,  it  was  the  object  (i.e.  pulse-shape)  that  was  constrained  making  the 
phase  retrieval  unnecessary.  Constraints  on  the  estimates  are  required  to  avoid  the  trivial 
solution  where  the  object  is  the  data  itself  and  the  PSF  is  a  delta  function.  Finally,  the  pixel 
bias  must  be  estimated.  In  order  to  estimate  the  pixel  bias,  the  second  set  of  complete  data, 
qk  (x,  y ),  is  utilized.  The  Q  function  for  the  pixel  bias  is 


k  x  Y 

k= 1  x=l  y=  1 


djk  (x,  y )  Bold  (x,  y ) 
illd  (x,  y)  +  Bold  (x,  y) 
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(5.59) 


Setting  dQs/dB  (xQ,  y0)  =  0  and  solving  for  the  pixel  bias  results  in 
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(5.60) 


GEM  iterations  continue  and  cease  when  the  mean-square  error  (MSE)  violates  the  stopping 
criteria.  Once  the  stopping  criteria  is  reached,  range  estimates  are  determined  by  using  the 
NCC  method  on  the  object  estimate. 

The  object  estimation  steps  are: 


1 .  Initialize  object,  PSF,  and  pixel  bias 

2.  Perform  one  GEM  iteration  using  Equations  (5.53),  (5.58),  and  (5.60) 

3.  Determine  MSE  and  compare  to  stopping  criteria 

4.  Repeat  Steps  2  and  3  until  stopping  criteria  reached 

5.  Use  NCC  to  find  new  range  estimates  with  Equation  (4.17) 

The  initial  estimates  in  step  one  are  the  same  as  GEMP  with  the  exception  that  the  object  is 
initialized  by  a  matrix  of  ones. 


5.2  Simulation 

In  order  to  verify  the  theory,  a  simulation  scenario  was  developed  whereby  targets  are 
interrogated  by  a  3D  FLASH  LADAR  defined  by  the  parameters  from  Table  5.1.  The  goal 
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Table  5.1:  3D  FLASH  LADAR  parameters 


Parameter 

Value 

Detector  Array 

128  x  128 

Aperture  Diameter  (D) 

2  mm 

Mean  Wavelength 

1.55  //m 

Focal  Length 

0.30  m 

Target  Range 

5.21  m 

Transmit  Energy 

10  mJ 

Pulse  Standard  Deviation  (aw) 

3  ns 

Beam  Divergence 

0.009  radians 

Detector  Spacing 

100  fim 

Detector  Array  Fill  Factor 

100% 

Detector  Bandwidth 

0.5  /im 

Target  Reflectivity 

10% 

Solar  Irradiance 

10  Watts/ni2///m 

D/ra  Seeing  Condition 

1.43 

Frame  Rate 

30  Hz 

Time  Samples 

20 

Sample  Period 

1.876  ns 

is  to  improve  range  estimation  given  the  noisy,  blurry  data  observations.  Results  show  range 
estimation  improvement  by  performing  object  recovery  either  via  a  Wiener  filter  method  or 
GEM  algorithms  as  outlined  in  Sections  5.1.2  and  5.1.3.  Previous  research  has  taken  the 
approach  to  use  a  Wiener  filter  on  each  individual  range  slice  and  then  use  a  pixel-based 
ranging  method  on  the  resulting  “deblurred”  data  cube  [55].  The  PSF  for  the  Wiener  filter 
is  set  as  the  diffraction-limited  PSF  of  the  system.  Performance  will  illustrate  that  the  GEM 
algorithms  provide  increased  error  performance  over  the  Wiener  filter  while,  at  the  same 
time,  being  more  robust.  Again,  the  GEM  algorithms  are  more  robust  in  that  they  do  not 
need  to  know  the  point  spread  function,  unlike  the  Wiener  filter  technique. 

Using  a  Gaussian  transmitted  pulse,  a  3D  FLASH  LADAR  imaging  scenario  is  devel¬ 
oped  in  simulation  using  various  geometrical  shapes  as  targets  shown  in  Figure  5.2(a)-(f). 
One  important  clarification  on  the  receiver  optics  is  that  the  detector  array  has  an  effective 
fill  factor  of  100%  by  placing  a  micro-lens  array  in  front  of  the  pixels  to  focus  the  light  onto 
the  pixel.  Also,  the  data  includes  effects  from  an  average  atmospheric  turbulence  to  enable 
blind  deconvolution.  Range  estimates  are  also  determined  without  processing  to  enable 
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pixel  pixel  pixel 


Figure  5.2:  True  ranges  for  simulation  targets:  (a)  three  bars,  (b)  Many  bars,  (c)  Various 
blocks,  (d)  Cylinder,  (e)  Slanted  boards,  and  (f)  Connected  blocks.  The  target 
names  in  this  caption  correspond  to  the  targets  in  Table  5.2.  The  three  bar  tar¬ 
get  is  also  the  experimental  data  target.  Other  targets  illustrate  the  robustness 
of  the  estimation  algorithms. 
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pixel 


further  comparison  between  no  processing  and  object  recovery  attempts.  Results  for  all  the 
targets  and  methods  with  error  metrics  are  summarized  in  Table  5.2.  The  numbers  in  bold 
indicate  the  best  performer  for  the  data  set. 
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Figure  5.3:  Estimated  ranges  for  simulation  targets:  (a)  No  processing  -  three  bars,  (b) 
GEM0  processing  -  three  bars,  (c)  No  processing  -  Many  bars,  and  (d)  GEM0 
processing  -  Many  bars.  Utilizing  the  GEM0  algorithm,  simulation  results 
show  the  image  quality  improvement  and  improved  range  estimation  (RMSE 
improves  75%  for  3  bar  target). 
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Table  5.2  clarifications:  “RMSE”  is  root  mean  square  error  (RMSE)  in  meters  be¬ 


tween  the  true  ranges  and  estimated  ranges  of  a  target  and  is  calculated  by 


RMSE 


\ 


M  N  ,  „  N  2 

E  (R(m,n)  —  R(m,n) 

771=1  71=1  ' 

~MN 


(5.61) 


where  R  (m,  n )  are  the  true  ranges  and  R  (m,  n )  (Equation  (4.17))  are  the  estimated  ranges. 
“Corr”  is  an  image  quality  metric  referring  to  the  correlation  coefficient  between  the  true 
range  image  and  estimated  range  image  signifying  linear  relationship  strength  (not  to  be 
confused  with  the  NCC  method)  and  mathematically  give  by 
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(5.62) 


where  hr  and  &r  are  the  mean  and  standard  deviation  of  the  true  range  image  respectively 
and  (1r  and  crr,  are  the  mean  and  standard  deviation  of  the  estimate  range  image  respec¬ 
tively.  “OD”  refers  to  the  original  data  (OD)  with  no  deblurring  and  ranges  estimated  by 
the  NCC  method,  “WF”  relates  to  range  estimation  using  a  Wiener  Filter  technique  with 
NCC  [55],  “GEMp”  is  the  pulse-shape  estimation  GEM  algorithm,  and  “GEM0”  is  the 
object  estimation  GEM  algorithm. 

The  targets  of  primary  interest  are  the  three  bar  target  and  the  multiple  bar  target  be¬ 
cause  the  three  bar  target  is  also  the  experimental  target  and  the  multiple  bar  target  is  most 
sensitive  to  spatial  blurring  of  all  the  targets.  The  bar  targets  are  constructed  in  simulation 
consisting  of  two  flat,  perpendicular  optically  rough  surfaces  at  different  ranges.  Referring 
to  Figures  5.2(a)  and  (b),  the  first  surface  in  range  has  rectangular  cut-out  shapes  while  the 
second  surface  contains  no  cutouts.  This  type  of  target  was  chosen  to  highlight  not  only  the 
coupling/blurring  effects  of  the  pixels  along  the  edges  of  the  rectangles,  but  also  the  decou¬ 
pling  and  ranging  capability  of  the  GEM  algorithm.  The  other  targets  are  built  in  similar 
manner.  Bar  target  shapes  were  used  because  the  distances  and  shape  dimensions  can  be 
physically  measured  in  a  laboratory  environment  to  show  range  estimation  improvement. 
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Figure  5.4:  (a)  Using  the  data  from  Figure  5.3(c)-(d),  investigating  pixel  (8,23)  shows 

the  estimated  waveform  (object  plus  pixel  bias)  closely  matching  the  true 
waveform  while  the  detected  waveform  does  not.  The  estimated  range  is  6.6 
m  while  the  true  range  is  6.7  m.  The  algorithm  also  implicitly  estimates  the 
pixel  bias  term  accurately. 

(b)  Again,  using  the  data  from  Figure  5.3(c)-(d),  investigating  pixel  (17,14) 
shows  the  estimated  waveform  improving  upon  the  detected  waveform,  but 
not  able  to  match  the  true  waveform  as  well  as  the  previous  pixel.  The  esti¬ 
mated  range  is  5.7  m  while  the  true  range  is  6.7  m.  Incorrect  range  estimation 
after  the  GEM0  algorithm  relates  to  blurring  severity  (edges  of  cut-outs  in 
first  surface)  and/or  a  particularly  noisy  realization  from  the  Poisson  distri¬ 
bution. 
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Table  5.2:  Range  estimation  results  for  simulation  data 


Data  set 

Metric 

OD 

WF 

GEMP 

GEM0 

Three  bars 

RMSE  (m) 
Corr 

0.402 

0.767 

0.346 

0.830 

0.163 

0.963 

0.100 

0.984 

Many  bars 

RMSE  (m) 
Corr 

0.596 

0.687 

0.561 

0.664 

0.346 

0.786 

0.365 

0.794 

Slanted  boards 

RMSE  (m) 
Corr 

0.225 

0.945 

0.171 

0.971 

0.161 

0.967 

0.131 

0.983 

Cylinder 

RMSE  (m) 
Corr 

0.184 

0.877 

0.153 

0.925 

0.160 

0.945 

0.153 

0.962 

Various  blocks 

RMSE  (m) 
Corr 

0.473 

0.595 

0.209 

0.931 

0.344 

0.725 

0.175 

0.955 

Connected  blocks 

RMSE  (m) 
Corr 

0.208 

0.853 

0.133 

0.955 

0.158 

0.918 

0.112 

0.970 

Table  5.2  and  the  range  images  from  Figures  5.2  and  5.3  show  the  negative  effects  of 
the  blurring  on  range  estimation  juxtaposed  with  the  positive  effects  from  attempting  to  re¬ 
cover  the  original  object  through  Wiener  filtering  or  the  GEM  algorithms.  Figure  5.4(a) 
shows  pixel  waveforms  successfully  recovered  while  Figure  5.4(b)  exhibits  a  situation 
where  the  recovery  was  not  as  successful.  Implicit  in  the  results  is  the  ability  to  accurately 
estimate  the  pixel  bias.  Without  it,  the  object  model  falls  apart  and  range  error  becomes 
extremely  large. 

An  additional  concern  is  assessing  the  computational  times  for  the  object  retrieval 
methods  (WF,  GEMp,  and  GEM0).  Although,  it  should  be  noted  again  that  the  WF  method 
requires  the  PSF  to  be  known  a  priori.  The  computational  times  were  analyzed  using  opera¬ 
tions  counting.  For  example,  this  counting  means  that  an  addition,  divide,  or  multiplication 
count  as  one  operation  each.  Also,  the  Fast  Fourier  Transform  (FFT)  is  utilized  to  accom¬ 
plish  convolution  and  correlation  and  counts  as  Nlog2  (N)  operations  where  N  is  the  number 
of  points  [13].  The  number  of  operations  required  in  the  WF  method  (Section  5.1.1)  are 
6  x  105.  Implementing  steps  2-7  from  Section  5.1.2  until  the  stopping  criteria  is  reached, 
the  GEMp  algorithm  uses  1.8  x  10'  operations  per  iteration.  Finally,  the  GEM0  algorithm 
has  a  computation  burden  of  1.9  x  106  operations  per  iteration  while  performing  Step  2-4 
(until  stopping  criteria  violation)  from  Section  5.1.3.  The  numbers  computed  for  all  the 
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methods  correspond  to  a  pixel  area  of  40  x  40  with  18  range  samples.  One  cube  was  used 
for  the  WF  and  GEMP  methods  while  two  cubes  was  used  for  GEM0.  The  additional  bur¬ 
den  on  GEMp  is  from  the  inner  iterations  which  include  100  GEM  iterations  followed  by 
an  update  on  the  pulse-shape  and  a  repeat  of  the  process.  While  the  increase  in  computation 
time  is  substantial  compared  to  the  WF  method,  the  GEMP  and  GEM0  algorithms  represent 
a  significant  increase  in  capability  with  respect  to  range  accuracy  and  to  required  a  priori 
information. 

Through  simulation,  the  model  and  object  recovery  attempts  have  been  verified.  The 
final  step  is  to  use  experimental  data  to  validate  simulation  results. 

5.3  Experimental  Results 

Using  the  pre-processed  experimental  data  described  in  Chapter  III,  Table  5.3  and 
Figure  5.5  illustrate  the  range  estimation  benefits  of  object  retrieval.  The  bold  numbers  in 
the  table  indicate  the  best  performing  method  for  the  data  set.  The  pulse-shape  and  object 
estimation  give  an  RMSE  improvement  of  25%  and  34%  respectively  over  the  original  data. 
Additionally,  the  pulse-shape  and  object  estimation  give  an  RMSE  improvement  of  7%  and 
18%  respectively  over  the  Wiener  filter  algorithm.  Figure  5.5(c)  shows  the  image  quality 
improvement  over  the  original  data  range  image  in  Figure  5.5(b).  Pixel  waveforms  provide 
additional  information  on  the  object  recovery  abilities.  Figure  5.5(d)  demonstrates  this 
ability  on  a  second  surface  pixel,  (32,18),  where  the  raw  waveform  results  in  an  incorrect 
range  determination.  In  contrast,  the  object  recovery  algorithm  (GEM0)  yields  an  improved 
range  estimate  by  sufficiently  estimating  the  true  waveform.  Additionally,  attempts  were 
also  made  to  use  asymmetrical  pulses  in  the  NCC  method.  However,  range  estimation 
error  did  not  significantly  change.  This  result  is  partly  due  to  the  gain  correction  where  the 
waveform  becomes  more  symmetrical  after  correction.  Since  the  exact  pulse  emitted  by 
the  laser  is  unknown,  the  pulse-shape  estimation  used  the  best  approximation  based  on  the 
observed  data. 
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pixel 


Figure  5.5:  Experimental  target :  (a)  True  ranges  with  first  surface  at  5.21  m  and  second 
surface  at  6.43  m  with  1.22  m  of  separation  in  between  surfaces, 

(b)  Ranges  using  NCC  without  using  object  recovery 

(c)  Estimated  ranges  using  GEM0  algorithm  followed  by  NCC 

(d)  Considering  pixel  (32,18),  its  estimated  waveform  (object  plus  pixel  bias) 
shows  similar  results  from  the  simulated  data.  The  estimated  waveform  more 
closely  resembles  the  true  waveform  with  the  range  close  to  range  sample  9. 
Also,  the  algorithm  correctly  estimates  the  pixel  bias  confirming  that  the  bias 
must  arise  from  a  noise  source  following  the  Poisson  distribution  (i.e.  dark 
current). 


Table  5.3:  Range  estimation  results  for  experimental  data 


Data  set 

Metric 

OD 

WF 

GEMP 

GEM0 

3  bars 

RMSE  (m) 
Corr 

0.301 

0.818 

0.243 

0.883 

0.226 

0.900 

0.198 

0.924 
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5.4  Conclusions 


Utilizing  waveform  sampling  capability,  the  positive  effects  of  object  recovery  in  3D 
FLASH  LADAR  range  estimation  is  clearly  evident.  The  innovative  3D  FLASH  LADAR 
sensor  provides  both  an  imaging  and  ranging  ability  enabling  established  theory  to  be  ap¬ 
plied  to  a  novel  manner.  Given  simulation  and  experimental  results,  it  is  clear  the  chosen 
model  and  noise  sources  are  an  appropriate  choice  for  3D  FLASH  LADAR  data  operating 
under  certain  conditions  (“sular  mode”  meeting  spatial  sampling  requirements).  The  raw 
data  coming  off  the  sensor  does  not  fit  the  model,  but  straight-forward  pre-processing  steps 
convert  the  data  to  an  acceptable  form  for  the  algorithms. 

In  mild  spatial  blurring  conditions,  simulation  results  predict  that  the  GEM  algo¬ 
rithms  increase  range  estimation  performance  substantially  over  no-processing  and  the 
Wiener  filter  method.  Again,  the  Wiener  filter  even  has  an  unfair  advantage  because  it 
is  provided  with  the  exact  (or  estimated)  PSF  function  used  in  generating  the  data  while  the 
GEM  algorithms  have  to  estimate  the  PSF  Considering  the  experimental  data,  its  perfor¬ 
mance  is  nearly  diffraction-limited  as  evidenced  by  the  experimental  PSF  and  OTF.  How¬ 
ever,  the  GEM  algorithms  still  increase  range  estimation  performance  over  the  Wiener  fil¬ 
ter.  Supported  by  simulation  results,  it  is  appropriate  to  say  that  the  GEM  algorithm  would 
show  even  better  range  estimation  performance  versus  the  Wiener  filter  in  severe  isoplanatic 
atmospheric  blurring  conditions  or  with  sub-optimal  optics. 

A  trade-off  exists  for  Wiener  filter  and  object  recovery  algorithms  between  compu¬ 
tation  cost  and  range  accuracy.  The  Wiener  filter  is  the  least  computationally  taxing  object 
recovery  algorithm,  but  is  the  least  accurate  and  requires  a  priori  knowledge  of  the  PSF.  The 
GEM  algorithms  are  computationally  expensive,  but  provide  the  best  range  performance 
and  can  perform  blind  deconvolution.  Considering  the  GEM  algorithms,  the  pulse-shape 
estimator  is  extremely  valuable  in  that  it  can  perform  range  estimation  on  a  single  cube 
thereby  removing  potential  for  any  registration  or  timing  errors.  If  multiple  cubes  are  avail¬ 
able  and  properly  registered,  object  estimation  is  undoubtedly  the  best  algorithm  to  use. 
Although,  experimentally,  none  of  the  algorithms  were  able  to  match  the  success  found 
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in  simulation.  Any  residual  error  in  the  experimental  results  can  be  attributed  to  system 
noise,  the  detected  light  containing  residual  laser  speckle,  residual  gain  error,  and  detector 
blurring. 

There  are  prospective  avenues  for  continued  investigation  and  improvement.  The 
pulse-shape  estimation  is  very  dependent  on  the  selected  waveform  model.  Improvements 
in  the  range  estimation  would  be  realized  if  a  true  waveform  model  for  the  transmitted 
laser  pulse  was  derived  or  calculated  experimentally.  Errors  in  the  experimental  data  re¬ 
sult  from  assuming  a  generalized  shape  that  is  corrupted  by  distorting  effects  (spatial  blur, 
pixel  blur,  and  noise).  In  addition,  the  variable  of  interest  (range  term)  would  ideally  be 
directly  estimated.  The  maximum  likelihood  solution  for  the  range  term  could  be  achieved 
if  another  model  was  discovered.  The  algorithm  in  this  paper  extracts  the  range  from  the 
maximum  likelihood  solution  for  the  pulse-shape.  Also,  even  after  the  pre-processing  steps, 
the  experimental  data  exhibits  noisy  behavior.  A  more  thorough  characterization  of  the  3D 
FLASH  LADAR  noise  sources  would  augment  or  verify  the  chosen  noise  sources.  Finally, 
isoplanatic  imaging  is  valid  for  the  experimental  set-up  in  the  laboratory.  However,  object 
recovery  from  3D  FLASH  LADAR  observations  subject  to  heavy  anisoplanatic  turbulence 
would  provide  an  ability  to  improve  range  estimation  in  a  variety  of  field  or  operational 
situations. 
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VI.  Range  Separation  Performance  and  Optimal  Pulse-Width 
Prediction  of  a  3D  FLASH  LADAR  using  the  Cramer- Rao  Bound 


The  purpose  of  this  chapter  is  to  use  established  estimator  variance  theory  applied  to 
a  three  dimensional  FLASH  LAser  Detection  And  Ranging  (3D  FLASH  LADAR) 
sensor  in  order  to  bound  estimator  performance  and  enable  prediction  of  optimal  resolu¬ 
tion  performance.  The  theory  in  this  chapter  supports  “hit  mode”  and  “sular  mode”  of  3D 
FLASH  LADAR  operation.  Additionally,  pixel  spatial  and  temporal  integration  is  a  com¬ 
mon  concern  for  photo-detectors.  For  this  paper,  the  pixel  spatial  response  is  assumed  to 
be  ideal  while  the  temporal  integration  is  assumed  to  be  less  than  half  of  the  range  interval. 
Supplemental  simulations  were  completed  that  investigated  integration  effects.  The  results 
show  that  these  effects  have  a  negligible  effect  on  the  blurring  severity  and  received  pulses. 
The  entire  chapter  contains  novel  and  original  efforts:  CRB  derivation  on  range  and  spatial 
separation  as  well  as  target  amplitudes,  optimal  pulse-width  determination  given  simple 
and  complex  targets,  and  optimal  pulse-width  determination  given  a  normalized  pulse  def¬ 
inition. 

The  Cramer-Rao  Bound  (CRB)  is  an  important  theoretical  result  in  estimation  theory 
that  can  be  applied  to  numerous  fields  in  science  and  engineering  [84].  Pertaining  to  a  two 
point  target  scene  illuminated  by  a  3D  FLASH  LADAR  (Section  4.4.1),  the  CRB  is  utilized 
to  bound  the  range  separation  estimation  variance.  The  simple  scene  is  adopted  to  allow 
for  closed-form  results  and  to  allow  conclusions  to  be  drawn  about  the  effects  of  range 
separation  on  the  bound.  Once  the  range  separation  CRB  is  derived,  an  unbiased  range 
separation  estimator  is  developed  to  enable  comparisons  to  the  CRB.  The  expected  results 
are  shown  in  an  example  comparing  the  estimator  variance  and  the  bound  across  possible 
range  separations. 

Additionally,  the  CRB  is  used  to  predict  system  performance  which  could  aid  in 
the  LADAR  development  process.  Per  conventional  RADAR  theory,  range  resolution  im¬ 
proves  (i.e.  becomes  smaller)  as  the  effective  pulse-width  is  shortened  [76].  Although,  the 
RADAR  engineer  must  be  concerned  about  other  factors  as  well  to  include  the  high  peak 
power  requirements  of  a  narrow  pulse.  In  the  case  of  3D  FLASH  LADAR,  there  is  the  abil- 
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ity  to  produce  ultra-short  laser  pulses  in  the  femtoseconds  (10-15)  compared  to  the  laser 
pulse  in  the  nanoseconds  (10~9)  used  in  this  research  [66],  [88].  Along  with  benefits  to 
target  ranging  and  identification,  one  would  expect  that  the  increase  in  the  range  resolution 
would  be  improved  by  several  orders  of  magnitude  with  an  ultra-short  laser  pulse.  How¬ 
ever,  similar  to  the  RADAR  engineer  concerns  about  high  peak  power  for  short  pulses,  the 
LADAR  engineer  has  to  be  concerned  that  the  receiver  electronics  can  sufficiently  sample 
the  returned  pulse.  With  the  laser  pulse-width  lasting  in  the  tens  of  nanoseconds,  the  current 
receiver  technology  can  only  generate  a  finite  amount  of  samples  due  to  the  complicated 
design  that  is  required  to  sample  the  pulse  every  couple  of  nanoseconds  [19].  Recognizing 
the  design  issues,  CRB  theory  is  employed  to  analyze  the  trade-off  between  laser  pulse- 
width  and  range  sampling  interval.  CRB  theory  and  subsequent  simulation  determine  that 
there  is  an  optimal  pulse-width  that  produces  an  optimal  range  resolution  for  a  particular 
range  sampling  capability. 

The  chapter  is  organized  as  follows:  Section  6. 1  derives  the  range  separation  CRB  for 
a  two  point  target  scene,  Section  6.2  discusses  the  results  of  the  unbiased  range  separation 
estimator  and  compares  it  to  the  CRB,  and  Section  6.3  uses  CRB  and  simulation  to  find 
an  optimal  pulse-width  for  several  different  range  sampling  scenarios  for  the  two  target 
case.  This  section  also  finds  an  optimal  pulse-width  considering  more  complex  targets. 
Additionally,  an  optimal  pulse-width  for  the  two-target  scene  if  found  using  a  normalized 
pulse  that  is  independent  of  range  sampling.  Finally,  Section  6.4  draws  conclusions  based 
on  the  observed  results. 

6.1  CRB  on  Range  Separation  Estimation 

In  this  section,  the  CRB  for  range  separation  A*,  is  derived  using  the  two-point  target 
data  model  from  Section  4.4.1.  Other  bounds  are  determined  as  well  including  spatial 
separation  and  the  target  amplitudes.  For  a  particular  imaging  scenario,  the  range  separation 
CRB  is  shown  in  a  figure  across  the  possible  range  separations. 

For  multiple  unknowns,  the  CRB  is  defined  by  the  diagonals  of  the  Fisher  Information 
Matrix  (FIM)  inverse  and  provide  a  lower  bound  on  the  variance  of  any  unbiased  estimator 
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which  is  shown  for  a  general,  multiple  unknown  parameter  case  by  [84] 


var 


(6.1) 


where  “var”  is  the  variance,  (1-,  (D)  is  the  estimate  of  a  particular  unknown,  non-random 
variable  6(,  D  is  the  observation  space,  and  J  is  the  FIM.  The  elements  of  the  FIM  are 
the  negative  expected  value  of  the  double  derivative  log-likelihood  function  and  provides  a 
measure  of  the  amount  of  information  of  an  unknown  parameter  contained  in  the  random 
process.  Mathematically,  the  FIM  is  defined  by  [84] 


Jij  —  — E 


d2  In  P  (dk  (x,  y )  =  Dk  {x,  y)  Wk,  x,  y ) 

86,89 j 


(6.2) 


where  E  is  the  expected  value  operation,  “In”  is  the  natural  log,  and  P  is  the  probability 
mass  function  (PMF)  for  all  3D  FLASH  LADAR  observations  with  dk  (x.  y )  as  the  real¬ 
izations  of  the  observations.  Assuming  statistical  independence  of  each  volume  element 
(voxel),  the  PMF  for  the  data  model  is  defined  by 
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(6.3) 


where  the  assumed  dominant  noise  source  is  photon  (shot)  noise  described  by  the  Poisson 
distribution.  While  lasers  exhibit  partial  coherence  meaning  the  negative  binomial  distri¬ 
bution  should  be  used  for  the  light  statistics,  this  photon  noise  assumption  is  valid  when 
the  operating  environment  produces  a  large  enough  speckle  parameter  so  that  the  nega¬ 
tive  binomial  distribution  approaches  the  Poisson  distribution  [24].  Previous  3D  FLASH 
LADAR  work  has  shown  the  speckle  parameter  to  be  adequate  to  assume  the  Poisson  dis¬ 
tribution  [9].  Additionally,  the  Poisson  distribution  CRB  provides  a  lower  bound  to  the 
negative  binomial  CRB  considering  the  higher  negative  binomial  variance  [60].  This  fact 
creates  a  true  lower  bound  (most  pessimistic)  with  the  Poisson  distribution  CRB  regardless 
of  the  imaging  conditions. 
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After  performing  the  required  operations,  the  general  solution  for  the  FIM  elements 


is  determined  to  be  [9],  [39] 
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Particular  to  this  work,  there  are  four  non-random  unknown  variables  in  the  data  model, 
9  =  [Am,  Ak,  At,  Ar\,  resulting  in  a  4x4  FIM  with  its  elements  determined  to  be 
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where 
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NOTE:  The  “x”  in  J12  is  a  multiply  operation.  The  FIM  is  inverted  and  the  CRB  for  each 
of  the  unknowns  is  on  the  diagonal  of  the  inverted  FIM  matrix  with  the  range  separation 
CRB  at  [J_1]22-  The  purpose  behind  supplying  the  FIM  element  expressions  was  to  provide 
enough  information  to  enable  the  work  to  be  reproduced.  Although  an  example  plot  is  given 
later  in  the  section,  the  range  separation  CRB  expression  itself  is  not  shown  due  to  its  length 
and  complexity. 

Besides  the  four  non-random  unknown  parameters,  the  CRB  also  depends  on  non- 
random  known  parameters  to  include  apt,  o>(,  and  ts.  In  order  to  view  a  useful  plot,  all 
other  unknown  and  known  factors  are  held  constant  while  the  Ak  is  stepped  from  the  be¬ 
ginning  to  the  end  of  the  range  extents.  Following  this  procedure,  Figure  6.1  shows  the 
range  separation  bound  for  a  specific  scenario  with  Am  =  1  pixel,  07,  =  3  pixels,  apt  =  3 
ns,  At  =  0.5  x  104  photons,  Ar  =  2  x  105  photons,  B  (x,y)  ~  A(750,38)  in  units  of 
photons,  and  range  sampling  ts  =  1.876  ns.  These  values  were  chosen  to  represent  a 
scenario  where  the  3D  FFASH  FADAR  interrogates  adjacent  targets  with  different  reflec¬ 
tivities  while  experiencing  significant  turbulence  in  the  atmosphere.  Furthermore,  the  bias 
definition  is  consistent  with  estimation  results  from  experimental  data. 
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Figure  6.1:  This  plot  shows  an  example  CRB  with  Am  =  1  pixel,  cx^  =  3  pixels, 

Cpt  =  3  ns,  At  =  0.5  x  104  photons,  and  Ar  =  2  x  105  photons.  The 
bound  behaves  appropriately  considering  the  variance  goes  up  as  the  sepa¬ 
ration  becomes  smaller  corresponding  to  the  notion  that  close-in  targets  are 
tougher  to  resolve.  The  peak  of  the  bound  occurs  when  the  range  and  spa¬ 
tial  coupling  are  at  their  maximum.  Further,  when  the  range  separation  near 
zero,  the  range  coupling  is  diminished,  but  the  bound  doesn’t  go  to  exactly 
zero  because  the  spatial  coupling  is  still  present. 
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The  shape  of  the  curve  in  Figure  6.1  reflects  the  negative  effects  of  the  range  and 
spatial  blurring  as  the  targets  become  closer.  Although,  the  effects  of  range  blurring  are 
minimized  when  the  targets  are  at  nearly  identical  ranges  and  the  bound  primarily  depends 
on  the  spatial  blurring.  Additionally,  the  increase  in  the  bound  past  ±2  meters  of  range 
separation  is  due  to  the  truncation  of  the  pulse  at  those  ranges.  An  assumption  in  the  bound 
derivation  is  a  fully  contained  pulse  within  the  range  extents.  The  impact  is  negligible 
considering  the  eventual  application  of  the  CRB  towards  range  resolution.  Targets  with 
±2  meters  of  range  separation  would  be  easily  resolved.  Changes  in  these  values  affect 
the  bound  in  a  predictive  manner.  For  example,  increasing  ah  and  apt  doesn’t  affect  the 
general  shape  of  the  range  separation  CRB,  but  it  does  increase  the  bound’s  magnitude  due 
to  increased  spatial  and  range  blurring  hampering  range  separation  estimation. 

More  specifically,  Figure  6.2  shows  several  examples  of  how  the  range  separation 
CRB  is  affected  by  changing  parameters  in  the  model  including  target  amplitude,  blurring 
severity,  and  spatial  separation.  Each  individual  figure  holds  all  other  parameters  constant 
and  plots  the  CRB  while  changing  one  parameter.  For  example,  Figure  6.2(a)  changes  the 
unknown  target  amplitude  At  while  keeping  all  other  parameters  constant.  Unless  otherwise 
noted,  the  standard  values  for  the  parameters  are  :  apt  =  3  ns,  5m  —  1  pixel,  ay,  =  3  pixels, 
ts  =  1.876  ns,  At  =  0.5  x  104  photons,  and  Ar  =  2  x  105  photons.  The  next  few  para¬ 
graphs  detail  how  the  changing  parameters  effect  the  range  separation  CRB.  The  parameter 
changes  affect  only  the  bound  values,  but  not  the  general  shape  of  the  bound. 

Figure  6.2(a)  -  At  effects.  As  the  unknown  target’s  amplitude  is  increased,  the  bound 
decreases  meaning  that  higher  SNR  values  of  the  unknown  target  aids  in  range  separation 
estimation  efforts.  (Inversely  proportional  to  bound) 

Figure  6.2(b)  -  Ar  effects.  Changing  the  known  target’s  amplitude  has  the  opposite 
effect  of  At.  As  the  Ar  amplitude  is  increased,  the  bound  also  increases  meaning  that 
range  separation  estimation  becomes  more  difficult  due  to  the  increased  blurring  between 
the  targets.  In  other  words,  estimating  the  range  separation  becomes  very  difficult  when 
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CRB(Afcj  _  o  CRB(Afc) 


Figure  6.2:  Effects  on  CRB(Afc)  when  changing  several  parameters  in  the  model  includ¬ 
ing  target  amplitude,  blurring  severity,  and  spatial  separation. 

(a)  At  -  inversely  proportional  to  bound 

(b)  Ar  -  proportional  to  bound 

(c)  ah  -  proportional  to  bound 

(d)  Am  -  inveresly  proportional  to  bound 
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considering  a  very  bright  known  target  next  to  a  much  dimmer  unknown  target  and  vice 
versa.  (Proportional  to  bound) 

Figure  6.2(c)  -  ah  effects.  As  the  blurring  severity  increases,  the  bound  also  increases 
meaning  that  more  blurring  (i.e.  more  coupling  between  the  targets)  hinders  range  separa¬ 
tion  estimation  performance.  (Inversely  proportional  to  bound) 

Figure  6.2(d)  -  Am  effects.  Finally,  as  the  spatial  separation  increases,  there  is  a  log¬ 
ical  corresponding  range  separation  estimation  performance  improvement  due  to  decrease 
in  coupling  between  the  targets.  (Inversely  proportional  to  bound) 

6.2  Range  Separation  Estimation  Results 

Using  the  model  governed  by  Equation  (6.3)  and  the  standard  parameter  values  from 
the  previous  section,  an  unbiased  range  separation  estimator  from  Section  4.4.2  is  applied 
to  enable  comparisons  of  the  range  separation  estimator  variance  to  the  CRB.  Other  pixel- 
based  range  estimators  are  available  including  peak  detection,  matched  filtering,  and  nor¬ 
malized  cross-correlation.  However,  in  this  two-target  scenario,  these  estimators  are  all 
biased  because  they  assume  that  there’s  only  one  target  per  pixel.  While  one  may  try  to 
deblur  the  data,  the  operation  will  not  be  totally  successful  and  some  bias  will  still  result. 
The  estimator  used  in  the  subsequent  sections  is  different  in  that  it  is  defined  as  having  two 
targets  per  pixel  thus  eliminating  the  bias. 

Prior  to  comparison,  the  estimator  must  first  be  determined  to  be  unbiased.  An  es¬ 
timator  of  an  unknown  parameter  is  unbiased  if  the  expected  value  of  the  estimator  is  the 
unknown  parameter  itself  (i.e.,  on  average,  one  expects  the  estimator  to  choose  the  true 
value  of  the  parameter  to  be  estimated)  [84].  In  terms  of  this  simulation,  the  estimator  is 
considered  to  be  unbiased  if  the  bias  squared  contribution  to  the  mean  square  error  (MSE) 
is  small  compared  to  the  range  variance  contribution.  This  relationship  results  from  the  fact 
that  MSE  equals  the  range  variance  plus  bias  squared.  Analytical  methods  to  determine  the 
bias  are  available,  but  graphical  nature  of  the  algorithm  prohibits  such  undertakings  requir- 
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ing  the  generation  of  statistics  based  on  many  instantiations  of  the  observed  data  through 
simulation. 

Therefore,  a  two  point  target  simulation  is  constructed  reflecting  the  observation 
model  defined  in  Equation  (6.3).  Shown  in  Figure  6.3,  the  simulation  results  include 
MSE,  bias  and  range  separation  variance.  As  expected,  Figure  6.3(c)  shows  the  bias  de¬ 
crease  as  the  iterations  increase  with  a  small  bias  left  after  the  last  iteration.  Referring  to 
Figures  6.3(a),  (b),  and  (c),  the  range  variance  dominates  the  MSE  and  the  two  point  target 
estimator  is  determined  to  be  unbiased. 

With  the  estimator  established  as  unbiased,  the  range  separation  variance  is  com¬ 
pared  to  the  CRB  to  observe  how  each  is  affected  given  changes  in  the  range  separation. 
Figures  6.1  and  6.3(b)  show  the  CRB  and  range  separation  variance  respectively.  Although, 
in  order  to  gain  more  insight  and  show  trends,  Figure  6.4  compares  the  CRB  and  range 
variance  in  the  same  plot  where  the  similar  behavior  is  now  evident.  In  fact,  the  estima¬ 
tor  range  variance  is  such  that  it  approaches  equality  with  the  bound.  This  equality  would 
make  the  estimator  efficient  [84] .  The  definition  of  an  efficient  estimator  relates  to  the  CRB . 
CRB  theory  states  that  any  unbiased  estimator  must  a  variance  equal  to  or  greater  than  the 
bound.  An  efficient  estimator  is  an  unbiased  estimator  whose  variance  equals  that  of  the 
bound.  Although  this  estimator  was  shown  to  be  unbiased,  it  is  not  theoretically  guaranteed 
to  be  efficient.  In  addition,  toward  the  edges  of  the  range  separation  ±  1.5  m,  the  bound 
should  theoretically  go  to  zero  like  the  variance  does,  but  it  doesn’t  because  the  Gaussian 
pulse  never  goes  to  zero.  This  non-zero  bound  can  be  ignored  since  estimation  is  easily 
performed  at  those  range  separations  and  can  also  be  mitigated  by  using  other  pulse  models 
such  as  a  negative  parabolic  that  equals  zero  until  the  pulse  is  received  [63]. 

6.3  Optimal  Pulse-width  Investigation 

Referring  to  Equation  (4.22),  crpt  controls  the  pulse-width  of  the  received  signal. 
Pulse  interactions  with  the  target  cannot  be  controlled,  but  the  transmitted  pulse- width  can 
be  factored  into  the  design  of  the  FADAR  system.  Figure  4.1(b)  shows  an  example  of  a 
pulse  with  apt  =  0.88  ns.  Following  standard  RADAR  theory,  a  smaller  effective  pulse- 
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Figure  6.3:  These  plots  show  the  range  separation  estimation  results  of  a  two  point  target 
data  model  simulation. 

(a)  Mean  square  error  (MSE)  between  the  truth  data  and  the  estimate. 

(b)  Range  separation  estimate  variance. 

(c)  Each  curve  is  a  bias  calcuation  for  a  different  Ak  over  many  trials.  At 
each  trial,  the  estimated  range  is  an  average  of  the  previous  estimated  ranges 
(i.e.  a  running  average). 

(d)  Bias  results  taken  from  the  last  trial  from  (c).  Comparing  (a),  (b),  and  (c), 
it  can  be  seen  that  estimator  is  unbiased  due  to  the  variance  dominating  the 
MSE. 
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Figure  6.4:  Taking  the  results  from  Figures  6.1  and  6.3(b),  this  plot  compares  the  CRB 
and  the  simulated  range  variance  showing  agreement  both  in  the  Cramer-Rao 
inequality  and  in  the  curve  shapes. 

width  is  desirable  due  to  its  ability  to  resolve  targets  closer  together.  However,  given  the 
discrete  nature  of  digital  sampling  employed  by  the  electronic  receivers  (denoted  by  ts  from 
Equation  (4.19)),  a  smaller  pulse- width  may  actually  degrade  performance  due  to  aliasing. 

Therefore,  the  CRB  is  used  to  predict  an  optimal  pulse-width  using  a  range  resolution 
metric.  Subsequently,  a  simulation  is  accomplished  to  validate  the  CRB  results.  Due  to  their 
separate  and  distinct  methods,  agreement  on  the  optimal  pulse-widths  between  the  CRB 
and  simulation  lends  confidence  to  the  results.  The  range-resolution  metric  is  defined  by 
comparing  the  square  root  of  the  CRB  (or  range  separation  variance  in  simulation)  with  the 
actual  range  separation,  Ak.  Referring  to  Figure  6.5,  the  location  where  the  values  equal  is 
defined  as  the  range-resolution  of  the  system.  In  other  words,  if  the  actual  range  separation 
is  within  one  standard  deviation  of  the  range  separation  estimate,  then  targets  would  not  be 
resolvable  due  to  the  estimation  uncertainty.  This  definition  implies  that,  on  average,  the 
estimator  would  be  able  to  resolve  targets  separated  further  while  targets  closer  together 
than  that  value  are  not  resolvable.  After  searching  over  many  pulse-widths,  the  pulse-width 
that  provides  the  best  range  resolution  (i.e.  lowest  value)  is  selected  as  optimal. 
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Figure  6.5:  This  figure  shows  an  example  plot  of  how  the  range  resolution  metric  is  de¬ 
termined.  The  circled  value  is  the  range  resolution  and  corresponds  to  the 
location  where  the  square  root  of  the  CRB  is  equal  to  the  range  separation. 
At  smaller  range  separations,  the  square  of  the  CRB  is  greater  than  the  sepa¬ 
ration  and  vice  versa. 


Table  6.1:  Optimal  pulse- width  results  -  two  point  target 


Range  Sampling  (ts) 

CRB,  opt  (ns) 

Simulation,  opt  (ns) 

0.6 tso  (1.126  ns) 

0.52 

0.52 

0.8fso  (1.500  ns) 

0.70 

0.70 

tso  (1.876  ns) 

0.88 

0.88 

1.2£so  (2.251  ns) 

1.06 

1.04 

1.4tso  (2.626  ns) 

1.22 

1.16 

6.3.1  Optimal  Pulse  for  Two  Point  Target.  Considering  the  two  point  target 
scenario,  Table  6.1  summarizes  the  optimal  pulse-width  apt  CRB  and  simulation  results 
for  several  range  sampling  cases  varying  from  faster  (0.6£so)  to  slower  (1.4 tso)  electronics. 
Figure  6.6  shows  the  data  points  for  each  range  sampling  case.  For  a  particular  ts,  changing 
the  pulse  to  be  either  narrower  or  wider  than  the  optimal  results  in  an  increase  in  the  bound 
or  variance  and  deterioration  in  the  range  resolution.  The  minimum  value  of  each  curve 
corresponds  to  the  reported  optimal  pulse- width. 
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apt  (ns) 
(a) 


Figure  6.6:  (a)  For  differing  range  sampling  cases,  the  range  resolution  derived  from  the 

CRB  is  plotted  versus  the  pulse-width.  As  the  range  sampling  ts  becomes 
either  faster  (0.6 tso  and  0.8 tso)  or  slower  (1.2tso  and  1.4iso),  the  optimal 
pulse-width  respectively  becomes  narrower  or  wider  with  a  corresponding 
improvement  or  degradation  in  the  range  resolution. 

(b)  The  simulation  range  resolution  determined  from  the  range  separation 
variance  is  plotted  versus  the  pulse-width.  As  expected,  the  resolution  values 
are  larger  than  those  predicted  by  CRB  theory.  Also,  the  optimal  pulse- width 
trends  in  a  similar  manner  as  the  CRB  results. 
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Figure  6.7:  Utilizing  the  optimal  apt  values  from  Table  6.1,  this  plot  shows  the  near 

exact  percentage  change  of  the  CRB  and  simulation  optimal  pulse-widths 
with  respect  to  the  percentage  change  in  range  sampling. 


As  can  be  seen  in  Figure  6.7,  the  optimal  pulse-width  scales  in  a  similar  manner  as 
the  range  sampling.  For  example,  if  the  range  sampling  was  reduced  by  80%,  then  the 
optimal  pulse- width  also  changed  by  approximately  80%  for  both  the  CRB  and  simulation 
results. 


6.3.2  Optimal  Pulse  for  Complex  Targets.  The  goal  of  this  section  is  to  show 
through  simulation  that  the  optimal  pulse- width  theory  holds  for  more  complex,  two  surface 
targets.  The  CRB  theory  in  Section  6.1  was  developed  for  a  simple,  two  point  target  and 
doesn’t  directly  pertain  to  these  new  targets.  However,  since  the  bound  and  the  ensuing 
simulation  both  predict  an  optimal  pulse-width  for  the  simple  two  point  target,  intuition 
dictates  that  an  optimal  pulse-width  could  also  be  found  for  more  complex  targets.  Three 
additional  targets  are  selected  (Section  5.2):  multi-bar,  three-bar,  and  connected  blocks. 
The  first  surface  of  all  these  targets  is  a  flat,  optically  reflective  board  with  shapes  cut  out. 
The  second  surface  is  also  flat  and  optically  reflective  and  is  placed  behind  the  first  board  at 
a  specific  range  separation.  This  second  surface  has  no  cut-outs.  Depending  on  the  target 
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Table  6.2:  Optimal  pulse-width  results  -  complex  targets 


Target 

Simulation,  apt  (ns) 

Multi-bar 

1.1 

Three-bar 

1.0 

Connected  blocks 

1.0 

geometry,  blurring  strength,  and  surface  reflectivity,  a  given  pixel  might  contain  significant 
contributions  from  either  one  surface  or  two  surfaces.  For  conciseness,  the  multi-bar  target 
method  and  results  are  discussed  in-depth  while  just  the  results  are  shown  for  the  other  two. 

The  method  to  determine  the  optimal  pulse-width  is  the  exactly  the  same  as  the  pre¬ 
vious  section:  vary  the  range  separations  and  accumulate  statistics  at  those  separations 
and  choose  the  pulse- width  that  produces  the  minimum  range  resolution.  Particular  to  this 
scenario,  the  pixel-based  two  surface  estimator  from  Section  4.5  is  used  to  generate  the 
estimates  using  a  threshold  of  7  =  0.97.  This  threshold  favors  the  “one  surface  pixel” 
model  due  to  false  peaks  created  by  noisy  realizations  of  the  incident  low-light  levels.  The 
first  surface  (in  range)  is  fixed  and  assumed  known  while  the  second  surface  is  placed  at 
successively  larger  distances  from  the  first  surface.  At  each  range  separation,  only  the  pix¬ 
els  classified  as  “two  surface”  are  used  in  order  to  keep  the  model  as  close  as  possible  to 
the  simple  two-point  target  CRB  of  Section  6.1.  The  estimation  statistics  collected  include 
variance,  mean  square  error,  and  bias.  Due  to  the  complexity  of  the  target  and  inherent 
coupling  between  adjacent  pixels,  low  light  levels  (15-30  received  photons)  where  required 
to  increase  the  effect  of  the  variance  on  the  observed  data. 

A  simulation  is  set  up  where  the  complex  targets  are  interrogated  by  a  3D  FLASH 
LADAR.  Results  of  the  simulation  are  shown  in  Table  6.2  and  Figures  6.8(b),  6.9(b), 
and  6.10(b)  where  the  optimal  pulse- width  standard  deviation  of  1.0  or  1.1  ns  show  moder¬ 
ate  agreement  to  the  CRB  results  (crpt  =  0.88  ns)  from  the  simple,  two  point  target.  Again, 
there  is  no  claim  that  the  results  have  to  match,  but  that  fact  that  they  are  close  for  several 
different  targets  is  encouraging. 

From  [5],  it  was  shown  that  the  two- surface  estimator  is  unbiased  given  a  simple 
scene.  However,  in  order  to  justifiably  compare  to  the  optimal  pulse-width  predicted  by 
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Figure  6.8: 


(a)  True  target  scene. 

(b)  Optimal  pulse  results  using  against  a  complex  target  with  ts  =  tso. 


Figure  6.9:  (a)  True  target  scene. 

(b)  This  plot  shows  the  optimal  pulse  results  for  a  three-bar  target  with  ts  = 

t  so- 
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Figure  6.10:  (a)  True  target  scene. 

(b)  The  optimal  pulse  is  shown  for  a  connected  blocks  target  with  ts  =  tso. 


the  CRB  in  Section  6.3.3,  the  estimator  must  be  shown  to  be  unbiased  given  the  complex 
scenes  where  the  convolution  effects  introduce  severe  bias  into  estimates.  Thus,  the  light 
levels  had  to  be  lowered  to  levels  where  the  maximum  peaks  of  the  observed  waveforms  are 
between  15  and  30  photons.  This  low  light  level  allows  for  variance  to  have  a  significant 
impact  on  observed  photon  counts. 

Figures  6.11(a)-(d)  and  6.12(a)-(d)  show  the  statistics  for  different  range  separations 
and  pulse-widths  respectively  for  the  multi-bar  target.  The  significant  factor  across  the 
plots  is  that  there  is  a  relatively  low  or  non-dominating  bias  in  the  region  of  range  separa¬ 
tion  where  the  range  resolution  is  selected.  This  region  for  best-performing  pulse-widths 
from  Figure  6.8(b)  ( apt  =  0.9, 1.0, 1.1,  and  1.2  ns)  shows  low  bias  and  therefore  variance 
dominance  in  the  MSE  near  the  selected  range  resolutions  of  0.8  to  0.9  centimeters  (re¬ 
ferring  to  Figure  6.12(a)-(d)).  These  results  indicate  that  some  regions  of  the  search  space 
(range  separation  and  pulse-width)  are  more  biased  than  others.  The  areas  that  are  of  most 
interest  tend  to  be  less  biased  and,  thus,  permit  a  conditional  comparison  to  the  CRB  results. 
Furthermore,  the  optimal  pulse-width  results  for  the  three-bar  and  connect  blocks  targets 
show  that  their  estimation  statistics  act  in  a  similar  manner.  With  the  variance  dominating 
in  the  areas  of  interest,  the  majority  of  the  error  lies  within  the  variance  allowing  for  the  use 
of  the  range  resolution  metric  and  subsequent  comparison  to  the  CRB  results. 
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Figure  6.11:  Statistics  across  pulse- widths  for  various  range  separations  for  the  multi-bar 
target.  The  statistics  include  mean  square  error,  bias,  and  range  variance. 

(a)  Range  separation  1.6  mm. 

(b)  Range  separation  11.6  mm. 

(c)  Range  separation  21.6  mm. 

(d)  Range  separation  31.6  mm. 
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Figure  6.12:  Statistics  across  range  separations  for  various  pulse-widths  for  the  multi-bar 
target.  The  statistics  include  mean  square  error,  bias,  and  range  variance. 

(a)  opt  =  0.9  ns. 

(b)  apt  =  1.0  ns. 

(c)  opt  =  1.1  ns. 

(d)  opt  =  1.2  ns. 
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6.3.3  Optimal  Pulse  for  a  Two  Point  Target  using  Normalized  Pulse  Definition  . 
Another  method  to  investigate  the  optimal  pulse- width  is  to  normalize  the  pulse  with  respect 
to  the  range  sampling.  Rather  than  define  the  width  of  the  pulse  by  crpt  in  units  of  seconds, 
the  pulse  is  defined  by  an  in  units  of  samples  given  by 

pW=7 ikexp{^r}  <67) 

where  n  €  [1,  Ar]  is  an  integer  range  sample  and  the  standard  deviation  is 


&pd  C&pt 

ct §  Ct  g 


®pt 

ts 


(6.8) 


with  opd  as  the  pulse  standard  deviation  in  units  of  meters,  c  is  the  speed  of  light  in  meters 
per  seconds,  ts  is  the  range  sampling  in  seconds  per  sample,  and  ap,  is  the  pulse  standard 
deviation  in  units  of  seconds.  The  benefit  of  this  definition  is  the  ability  to  determine  an 
optimal  pulse-width  independent  of  range  sampling  capability. 

Pertaining  to  the  range  resolution  metric,  the  optimal  pulse  standard  deviation  on  is 
found  by  using  the  same  investigation  procedures  as  the  previous  section.  Figure  6.13(a) 
shows  the  optimal  pulse  shape  in  terms  of  a  standard  deviation  measured  in  samples.  Using 
Equation  (6.8),  the  optimal  standard  deviation  in  seconds  opt  can  be  found  for  a  particular 
ts  as  seen  in  Figure  6.13(b).  These  values  for  opt  are  consistent  with  the  values  from  the 
previous  section. 

Considering  all  the  above  optimal  pulse-width  studies,  an  important  observation  is 
the  number  of  significant  samples  across  the  pulse  for  the  optimal  pulse-widths  from  Ta¬ 
bles  6.1,  Table  6.2,  and  the  normalized  pulse  section.  In  each  case,  the  number  of  significant 
samples  across  the  optimal  pulse  is  three.  Referring  to  Figure  4.1(b),  a  significant  sample 
is  defined  as  a  non-zero,  sizable  contributor  (>  2%  of  pulse  peak  value)  to  the  waveform. 
This  consistent  number  of  significant  samples  indicates  that  more  samples  (i.e.  larger  opt ) 
than  the  optimal  harms  the  range  resolution  capability  while  fewer  samples  fails  to  provide 
the  estimator  with  enough  data  and  under-samples  the  received  pulse  causing  the  variance 
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(a)  Using  the  normalized  pulse  model,  this  graph  shows  the  CRB  optimal 
pulse  standard  deviation  referring  to  Equation  (6.7). 

(b)  CRB  optimal  pulse  versus  range  sampling.  The  optimal  pulse  width 
changes  proportionally  as  the  sampling  changes. 
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to  increase.  Even  for  a  complex  target  whose  estimation  variance  is  not  described  by  the 
CRB  theory,  the  simulated  optimal  pulse-width  still  maintains  the  three  significant  samples. 

6.4  Conclusions 

The  CRB  is  used  to  bound  system  performance  by  giving  a  measure  of  how  well  an 
unbiased  estimator  can  perform.  In  this  paper,  the  CRB  bounds  the  performance  on  the 
ability  to  estimate  the  range  separation  of  two  spatial  point  targets.  Using  the  point  target 
estimator,  the  unbiased  estimator  variance  is  shown  to  be  bounded  by  the  CRB  with  both 
trending  in  a  similar  manner. 

Additionally,  the  CRB  can  be  used  to  not  only  predict  performance,  but  to  give  a 
LADAR  engineer  the  capability  to  predict  critical  LADAR  performance  without  regard 
to  the  particulars  of  estimation.  Considering  the  optimal  pulse-width  study,  all  the  range 
sampling  cases  produced  an  optimal  pulse- width  with  a  similar  number  of  significant  sam¬ 
ples  across  the  pulse.  Three  significant  samples  across  the  pulse  provides  the  best  range 
resolution  while  not  experiencing  the  ill-effects  of  under-sampling. 

The  overriding  conclusion  is  that  a  shorter  pulse-width  in  the  femtoseconds  does 
not  always  provide  improved  range  resolution  performance.  There  is  a  range  resolution 
performance  link  between  the  electrical  circuitry  sampling  capability  and  the  width  of  the 
transmitted  pulse.  In  conjunction  with  other  performance  goals  and  design  limitations,  the 
LADAR  engineer  concerned  about  range  resolution  should  not  just  focus  on  a  shorter  pulse- 
width  without  making  improvements  in  the  receiver’s  capability  to  sample  the  detected 
return  pulses. 
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VII.  Summary 


The  contributions  of  this  research  increase  the  body  of  knowledge  and  the  ultimate  per¬ 
formance  of  three  dimensional  FLASH  LAser  Detection  And  Ranging  (3D  FLASH 
LADAR).  By  decreasing  ranging  error,  characterizing  estimation  performance,  or  optimiz¬ 
ing  system  parameters,  sensor  capability  has  been  enhanced.  This  increase  is  aided  by  an 
increase  in  the  data  model  fidelity  from  previous  contributions  in  the  field  which  have  all  but 
ignored  the  spatial  blurring  effect  of  the  image  formation  process.  By  incorporating  these 
blurring  effects,  new  problems  and  opportunities  to  exploit  the  data  were  formed.  Solutions 
were  subsequently  found  that  took  advantage  of  the  “extra”  information  available. 

Given  3D  FLASH  LADAR’s  future  in  imaging,  computer  vision,  guidance,  naviga¬ 
tion,  and  targeting,  this  work  builds  a  basis  of  understanding  concerning  data  models  and 
data  processing.  The  research  has  unique  attributes  that  are  rare  considering  other  other  sig¬ 
nal  processing  research.  The  3D  FLASH  LADAR  is  a  powerful  sensor  combining  RADAR 
principles,  laser  transmission,  waveform  processing,  and  electro-optical  phenomena  and 
requires  equally  powerful  algorithms  to  estimate  and  characterize  parameters  of  interest. 
The  focus  of  this  research  centers  around  a  particular  parameter:  target  range.  Although, 
other  areas  of  research  may  be  possible  given  that  it  is  shown  that  a  particular  data  model 
appropriately  characterizes  data  from  a  3D  FLASH  LADAR  sensor. 

This  chapter  summarizes  the  previous  chapters  in  the  document,  reviews  the  signifi¬ 
cant  research  contributions,  and  outlines  several  avenues  for  future  research. 

7.1  Chapter  Summary 

Chapter  II  provided  background  theory  related  to  imaging,  coherence,  deconvolu¬ 
tion,  maximum  likelihood  principles,  and  generalized  expectation  maximization.  It  also 
presented  the  data  model  used  prevalently  in  this  dissertation.  Finally,  previous  research 
was  summarized  in  the  areas  of  3D  LADAR  hardware  and  data  processing,  blind  deconvo¬ 
lution,  performance  bounding,  and  parameter  optimization. 

Chapter  III  detailed  the  3D  FLASH  LADAR  hardware  and  laboratory  settings  used  in 
experimental  collects.  This  chapter  also  describes  the  procedures  used  to  condition  the  data 
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for  appropriate  use  for  the  selected  mathematical  model.  Finally,  it  explains  how  the  ex¬ 
perimental  point-spread-function  is  determined  as  well  as  the  speckle  parameter  estimation 
which  increases  confidence  with  the  chosen  incoherent  light  model. 

Chapter  IV  contained  the  pertinent  range  estimation  algorithms  to  include  peak  de¬ 
tection,  maximum  likelihood,  normalized  cross-correlation  (matched  filtering),  a  two  point 
target  range  and  spatial  separation  estimator,  and  a  two  surface  range  separation  estima¬ 
tor.  The  normalized  cross-correlation  is  usually  the  method  of  choice  due  to  its  ability  to 
perform  intersample  ranging  and  handle  truncated  waveforms.  Although,  it  may  encounter 
some  bias  depending  on  spatial  blurring  severity  and/or  the  effectiveness  of  deblurring  (de- 
convolution)  algorithms. 

Chapter  V  implemented  object  deconvolution  and  blind  deconvolution  on  blurry, 
noisy  data  observations  using  simulation  and  experimental  data  showing  that  object  re¬ 
covery  improves  range  estimation. 

Chapter  VI  derived  the  CRB  for  range  and  spatial  separation  as  well  as  amplitude 
estimation  considering  a  two-point  target  scene.  Considering  several  range  sampling  cases, 
the  range  separation  CRB  and  estimator  also  predicted  an  optimal  pulse-width  that  pro¬ 
vides  the  best  range  resolution.  Additionally,  an  optimal  pulse-width  is  found  using  more 
complex  targets  and  for  a  normalized  pulse-definition. 

7.2  Summary  of  Contributions 

7.2.1  Improving  Range  Estimation  by  Spatial  Processing.  The  benefits  of  spatial 
processing  are  evident  when  comparing  range  estimation  before  and  after  the  object  recov¬ 
ery  algorithms.  Simulation  shows  substantial  increase  in  the  image  quality  and  decrease 
in  the  range  estimation  errors.  The  experimental  data  shows  improvement,  but  not  as  dra¬ 
matic  due  primarily  to  the  excess  speckle  noise  evident  in  the  data.  A  favorable  result  is  that 
the  blind  deconvolution  methods  outperformed  deconvolution  even  when  the  deconvolution 
was  given  the  exact  form  of  the  blurring  function.  The  ability  to  process  three-dimensional 
data  in  two  dimensions  and  range  in  the  third  dimension  was  the  innovative  vision.  This 
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contribution  proves  that  this  vision  can  be  theoretically  proven  and  demonstrated  through 
simulation  and  experimentation. 

7.2.2  Unbiased  Two  Point  Target  Temporal  and  Spatial  Estimator.  Using  a  two- 
point  target  scene,  an  unbiased  spatial-temporal  estimator  is  derived  that  is  able  to  accu¬ 
rately  estimate  spatial  and  range  separations  and  received  amplitudes.  Through  simulation, 
its  range  separation  estimator  variance  compares  favorably  to  that  predicted  by  the  CRB. 
This  estimator  can  even  handle  truncated  waveforms  which  is  not  predicted  or  suitable  han¬ 
dled  by  the  bound.  Without  this  unbiased  estimator,  the  CRB  results  would  have  nothing 
to  compare  to  and  there  would  be  less  confidence  its  conclusions.  The  agreement  between 
the  CRB  and  simulation  variance  is  exceptionally  significant  and  vital  since  they  arrive  at 
virtually  the  same  answer  by  different  methods. 

7.2.3  Lower  Bound  on  Range  Separation  Estimation.  Through  the  use  of  a  two- 
point  target  scene,  the  CRB  on  range  separation  estimation  is  derived.  The  CRB  on  spatial 
separation  and  amplitude  estimation  is  derived  as  well.  Pertaining  to  range  separation,  the 
CRB  shows  that  range  separation  does  not  severely  affect  estimation  performance  until  the 
targets  are  close.  When  the  range  separation  is  identically  zero,  the  bound  does  not  go  to 
zero  due  to  the  spatial  blurring  still  present  in  the  data  (targets  are  still  spatially  separated). 
Additionally,  the  shape  of  the  bound  is  remains  constant  independent  of  several  parameters 
including  spatial  blurring  strength,  signal-to-noise  ratio  (SNR)  of  the  reflectors,  and  spatial 
separation.  The  dependence  between  the  bound  and  the  parameters  is  the  absolute  level  of 
the  CRB.  As  the  blurring  strength  increases,  the  bound  also  increases  and  vice  versa.  When 
the  SNR  of  the  unknown  target  reflector  is  either  increased  or  decreased,  the  bound  acts 
oppositely  and  decreases  or  increases  respectively.  These  results  should  be  intuitive  where 
increased  blurring  would  make  estimation  more  difficult  (i.e.  variance  would  increase)  and 
an  increase  in  the  unknown  target  SNR  would  enhance  the  estimation  abilities  (i.e.  variance 
would  decrease). 
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7.2.4  Optimal  Pulse-Width  based  on  Range  Resolution.  With  this  contribution, 
it  is  shown  that  the  CRB  can  perform  system  parameter  optimization  with  respect  to  a  very 
important  system  characteristic  and  requirement  -  range  resolution.  Independent  of  estima¬ 
tor  choice,  the  bound  shows  that  an  optimal  pulse- width  exists  whereby  the  expected  range 
resolution  is  minimized.  After  developing  an  unbiased  estimator  for  the  target  scene,  the 
optimal  pulse-widths  predicted  by  the  CRB  are  verified  through  simulation.  The  agreement 
between  CRB  and  simulation  is  very  significant  given  that  they  arrive  at  range  separation 
variance  either  through  Fisher  information  theory  or  through  repeated  trials  using  a  simu¬ 
lation. 

Furthermore,  the  range  sampling  interval  is  both  increased  (slower  electronics)  and 
decreased  (faster  electronics)  which  shows  the  resultant  optimal  pulse-shape  becoming 
wider  and  narrower  respectively.  In  other  words,  faster  electronics  that  sample  the  range  di¬ 
mension  faster  can  incorporate  a  narrower  pulse-width  and  achieve  better  range  resolution. 
To  lend  confidence  to  the  results,  optimal  pulse-widths  are  also  found  for  more  complex 
targets.  Also,  in  terms  of  samples,  an  optimal  pulse- width  using  the  CRB  is  found  using  a 
normalized  pulse  model.  This  definition  means  that  the  results  are  independent  of  the  range 
sampling  interval.  Finally,  and  perhaps  most  enlightening,  all  the  optimal  pulse-width  re¬ 
sults  reflect  that  the  received  pulse  needs  to  have  three  significant  samples  in  the  received 
data.  Fewer  significant  samples  caused  by  a  narrower  transmitted  pulse  or  target  interac¬ 
tions  does  not  provide  enough  information  to  match  the  model  pulse-shape  and  could  even 
be  entirely  missed  by  the  electronics.  Following,  while  providing  enough  information, 
more  significant  samples  would  certainly  be  less  optimal  by  degrading  range  resolution 
performance. 

7.3  Future  Research 

There  are  numerous  additional  research  avenues  available  with  respect  to  3D  FLASH 
LADAR  and  data  processing  including  the  following: 
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7.3.1  Outlier  Detection.  To  overcome  speckle  produced  by  coherent  lasers  and 
increase  signal-to-noise  ratio  (SNR)  because  of  low  return  signal  levels  due  to  low  trans¬ 
mit  energy,  propagation  distances,  or  target  reflectance,  3D  LADAR  systems  may  need  to 
average  many  cubes  of  scene  data.  In  practice,  certain  cubes  may  be  warped  due  to  mis¬ 
registration  or  atmospheric  effects.  Also,  a  particular  pixel  may  be  defective  for  a  small 
amount  of  time  resulting  in  an  out-of-family  pixel  that  is  “warped”  in  a  cube  of  data.  If 
the  warping  is  severe  enough,  the  averaging  process  may  be  negatively  affected  by  these 
particular  cubes  or  pixels.  It  would  be  advantageous  to  system  performance  (i.e.  object 
recovery,  range  estimation)  to  develop  a  means  to  intelligently  remove  these  frames  before 
the  averaging  process. 

7.3.2  FOliage  PENetration  (FOPEN)  Capability  Investigation.  A  key  military 
mission  for  any  imaging  or  ranging  sensor  is  the  ability  to  recognize  man-made  targets  un¬ 
der  foliage  that  can  either  be  man-made  itself  (camouflage)  or  natural  (trees).  Successful, 
experimental  efforts  have  already  been  accomplished  trying  to  ascertain  the  3D  FLASH 
LADAR’s  FOPEN  capability.  However,  a  rigorous  theoretical  model  has  not  been  adopted 
yet.  This  model  and  subsequent  simulation  and  experimental  investigation  would  numeri¬ 
cally  characterize  FOPEN  potential  in  a  myriad  of  environments  including  different  cam¬ 
ouflage  configurations,  look  angles,  weather  conditions  (e.g.  wind  velocity),  and  targets. 

7.3.3  Pixel  Impulse  Response  Deconvolution.  As  with  any  high-performance 
military  hardware,  characterization  under  environmental  operational  conditions  is  a  manda¬ 
tory  exercise.  The  operator  must  know  the  limits  where  one  would  expect  nominal  perfor¬ 
mance.  As  part  of  the  hardware  characterization  efforts,  the  pixel  impulse  response  impact 
on  the  reflected  pulse  is  important  when  developing  accurate  pulse  models.  The  pixel’s  im¬ 
pulse  response  is  not  ideal  and  does  have  some  distortion  effect  on  the  returned  waveforms. 
Using  simulated  and  experimental  data,  the  research  effort  would  calculate  the  distortion 
severity  and  dependence  on  system  parameters. 
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7.3.4  Object  blind  deconvolution  using  partially  coherent  light  model.  Described 
by  the  Poisson  distribution,  the  incoherent  light  model  used  in  this  research  is  an  approxi¬ 
mation  for  the  partially  coherent  model  which  is  the  most  accurate  portrayal  of  the  detected 
laser  light.  Detected  partially  coherent  light  is  statistically  described  by  the  negative  bino¬ 
mial  (NB)  distribution.  Blind  deconvolution  methods  using  the  NB  distribution  are  cum¬ 
bersome  and  nearly  intractable  unless  the  point- spread-function  can  be  parameterized.  If  an 
object  estimation  method  can  be  found  to  use  the  Generalized  Expectation  Maximization 
(GEM)  algorithm  with  the  NB  distribution,  the  resulting  estimator  would  theoretically  out¬ 
perform  the  object  estimator  in  this  research  due  to  the  increase  in  noise  model  accuracy. 
The  key  issue  in  the  derivation  comes  when  taking  the  conditional  expectation  of  the  log  of 
the  complete  data  with  respect  to  the  incomplete  data  and  the  old  estimates.  A  vital  prop¬ 
erty  of  the  Poisson  distribution  is  that  a  sum  of  Poisson  random  variables  is  still  Poisson. 
The  same  cannot  be  said  for  a  sum  of  NB  random  variables.  Consequently,  a  variation  of 
the  GEM  or  data  model  is  necessary  to  complete  the  derivation. 
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